双三次插值算法的C++实现与SSE指令优化

在上篇文章中，我们讲解了常见的最邻近插值算法、双线性插值算法和双三次插值算法的原理与实现，三种插值算法中双三次插值算法的插值效果最好，但其也是三种算法中计算复杂度最高、耗时最长的算法。本文在给出双三次插值C++代码的基础上，着重讲解如何使用SSE指令来优化该算法，并使用双三次插值来实现图像的缩放，比较SSE指令优化前后的耗时。

常见图像插值算法的原理与C++实现

1. 基于C++与Opencv的代码实现

算法原理在上篇文章中已经讲了，此处直接贴出代码：

float cubic_w_f(float x, float a)
{
  if (x <= 1)
  {
    return 1 - (a + 3)*x*x + (a + 2)*x*x*x;
  }
  else if (x < 2)
  {
    return -4 * a + 8 * a*x - 5 * a*x*x + a*x*x*x;
  }
  return 0.0;
}


void cal_cubic_coeff(float x, float y, float *coeff)
{
    float u = x - floor(x);
    float v = y - floor(y);


    u += 1;
    v += 1;
    float a = -0.15;


    float A[4];
    A[0] = cubic_w_f(abs(u), a);
    A[1] = cubic_w_f(abs(u - 1), a);
    A[2] = cubic_w_f(abs(u - 2), a);
    A[3] = cubic_w_f(abs(u - 3), a);


    for (int s = 0; s < 4; s++)
    {
      float C = cubic_w_f(abs(v - s), a);
      coeff[s * 4] = A[0] * C;
      coeff[s * 4 + 1] = A[1] * C;
      coeff[s * 4 + 2] = A[2] * C;
      coeff[s * 4 + 3] = A[3] * C;
    }
?}
  
uchar?cubic_inner(Mat?src,?float?x_float,?float?y_float,?float?a)
{
    float coeff[16];
    cal_cubic_coeff(x_float, y_float, coeff);  //计算权重系数


    float sum = 0.0;
    int x0 = floor(x_float) - 1;
    int y0 = floor(y_float) - 1;


    for (int i = 0; i < 4; i++)
    {
      for (int j = 0; j < 4; j++)
      {
        sum += coeff[i * 4 + j] * src.ptr(y0 + i)[x0 + j];
      }
????}
    return ((uchar)sum);
?}

2. SSE指令优化算法

首先，我们来看一下浮点型坐标点周围的4*4个整型点分别在x方向与y方向上与该浮点型坐标点的像素距离，假设浮点型坐标点的x坐标的小数部分为u，y坐标的小数部分为v，那么x方向与y方向上的距离如下图所示（每一格的像素距离为1）。

从左到右，x方向距离分别为1+u、u、1-u、2-u：

从上到下，y方向距离分别为1+v、v、1-v、2-v：

从而得到各个距离的取值范围：

1 ≤ dx0=1+u ≤ 2

0 ≤ dx1=u ≤ 1

0 ≤ dx2=1-u ≤ 1

1 ≤ dx3=2-u ≤ 2

1 ≤ dy0=1+v ≤ 2

0 ≤ dy1=v ≤ 1

0 ≤ dy2=1-v ≤ 1

1 ≤ dy3=2-v ≤ 2

双三次插值算法的权重计算公式为：

我们可以根据取值范围提前确定dxi与dyj的权重函数表达式（之前是分段函数），便于SSE指令的并行计算：

对于dx0、dx3、dy0、dy3，其权重函数表达式为：

对于dx1、dx2、dy1、dy2，其权重函数表达式为：

因此dx0、dx3、dy0、dy3的权重可以并行计算，dx1、dx2、dy1、dy2的权重同样也可以并行计算，假设浮点型坐标为(x, y)，权重的SSE指令并行计算代码如下：

float?u?=?x_float?-?floor(x_float);???//计算x坐标额小数部分
float?v?=?y_float?-?floor(y_float);???//计算y坐标额小数部分
float?a_mul_4?=?(a?+?a)?+?(a?+?a);???//提前计算权重公式中的4a
float?a_mul_5?=?a_mul_4?+?a;?????????//提前计算权重公式中的5a
float?a_mul_8?=?a_mul_4?+?a_mul_4;???//提前计算权重公式中的8a
float?a_add_3?=?a?+?3;   //提前计算权重公式中的a+3
float?a_add_2?=?a?+?2;   //提前计算权重公式中的a+2
__m128?a_m?=?_mm_set1_ps(a);??//a?a?a a
__m128?m_1?=?_mm_set1_ps(1.0);?//1.0?1.0?1.0?1.0
__m128?a_mul_4_m?=?_mm_set1_ps(a_mul_4);  //4a 4a 4a 4a 
__m128?a_mul_5_m?=?_mm_set1_ps(a_mul_5);??//5a?5a?5a 5a
__m128?a_mul_8_m?=?_mm_set1_ps(a_mul_8);  //8a 8a 8a 8a 
__m128?a_add_3_m?=?_mm_set1_ps(a_add_3);  //a+3 a+3 a+3 a+3
__m128?a_add_2_m?=?_mm_set1_ps(a_add_2);  //a+2 a+2 a+2 a+2


__m128?C30_A30?=?_mm_set_ps(2?-?v,?1?+?v,?2?-?u,?1?+?u);???//dy3?dy0?dx3?dx0
__m128?C21_A21?=?_mm_set_ps(1?-?v,?v,?1?-?u,?u);???//dy2?dy1?dx2?dx1


__m128?tmp0?=?_mm_sub_ps(_mm_mul_ps(a_m,?C30_A30),?a_mul_5_m);???//a*d?-?5a
tmp0?=?_mm_add_ps(a_mul_8_m,?_mm_mul_ps(C30_A30,?tmp0));???????//8a?+?d*(a*d-?5a)
tmp0?=?_mm_sub_ps(_mm_mul_ps(C30_A30,?tmp0),?a_mul_4_m);????//d*(8a?+?d*(a*d-?5a))?-?4a?=?w(dy3)?w(dy0)?w(dx3)?w(dx0)


__m128?tmp1?=?_mm_sub_ps(_mm_mul_ps(a_add_2_m,?C21_A21),?a_add_3_m);???//(a+2)*d?-?(a+3)
tmp1?=?_mm_mul_ps(_mm_mul_ps(C21_A21,?C21_A21),?tmp1);????//d*d*((a+2)*d - (a+3))
tmp1?=?_mm_add_ps(m_1,?tmp1);?????//1?+?d*d*((a+2)*d?-?(a+3))?=?w(dy2)?w(dy1)?w(dx2)?w(dx1)

以上代码运行之后得到权重如下（高位-->低位）：

tmp0：w(dy3) w(dy0) w(dx3) w(dx0)

tmp1：w(dy2) w(dy1) w(dx2) w(dx1)

全部的w(dxi)与w(dyj)都已计算完毕，但以上并不是我们想要的排列顺序，我们想要的排列顺序如下：

w(dy3) w(dy2) w(dy1) w(dy0)

w(dx3) w(dx2) w(dx1) w(dx0)

因此我们需要对tmp0与tmp1进行重新打包与排列：

__m128 A_m = _mm_unpacklo_ps(tmp0, tmp1);    //交替打包tmp0与tmp1的低位数据：tmp1[1] tmp0[1] tmp1[0] tmp0[0] = w(dx2) w(dx3) w(dx1) w(dx0)
__m128?C_m?=?_mm_unpackhi_ps(tmp0,?tmp1);????//交替打包tmp0与tmp1的高位数据：tmp1[3] tmp0[3] tmp1[2] tmp0[2]?= w(dy2) w(dy3) w(dy1) w(dy0)
A_m?=?_mm_shuffle_ps(A_m,?A_m,?_MM_SHUFFLE(2,?3,?1,?0));???//重新排列A_m中数据的顺序：w(dx3) w(dx2) w(dx1) w(dx0)
C_m?=?_mm_shuffle_ps(C_m,?C_m,?_MM_SHUFFLE(2,?3,?1,?0));???//重新排列C_m中数据的顺序：w(dy3) w(dy2) w(dy1) w(dy0)

接下来就可以计算W(i, j)=w(dxi)*w(dyj)了，由于i和j都取0、1、2、3，因此有4*4=16个W(i, j)，对应周围的4*4个整型点。代码如下：

__declspec(align(16)) float C[4];
?_mm_store_ps(C,?C_m);  //w(dy3) w(dy2) w(dy1) w(dy0)


?__m128?m128_C?=?_mm_set1_ps(C[0]);???//w(dy0) w(dy0) w(dy0) w(dy0)
?__m128?coeff0?=?_mm_mul_ps(A_m,?m128_C);  //W(3,0) W(2,0) W(1,0) W(0,0) = w(dx3)*w(dy0) w(dx2)*w(dy0) w(dx1)*w(dy0) w(dx0)*w(dy0)


?m128_C?=?_mm_set1_ps(C[1]);  //w(dy1) w(dy1) w(dy1) w(dy1)
?__m128?coeff1?=?_mm_mul_ps(A_m,?m128_C);  //w(dx3)*w(dy1) w(dx2)*w(dy1) w(dx1)*w(dy1) w(dx0)*w(dy1)


?m128_C?=?_mm_set1_ps(C[2]);   //w(dy2) w(dy2) w(dy2) w(dy2)
?__m128?coeff2?=?_mm_mul_ps(A_m,?m128_C);  //w(dx3)*w(dy2) w(dx2)*w(dy2) w(dx1)*w(dy2) w(dx0)*w(dy2)


?m128_C?=?_mm_set1_ps(C[3]);   //w(dy3) w(dy3) w(dy3) w(dy3)
?__m128?coeff3?=?_mm_mul_ps(A_m,?m128_C);  //w(dx3)*w(dy3) w(dx2)*w(dy3) w(dx1)*w(dy3) w(dx0)*w(dy3)

最后，就可以计算4*4个整型点的像素加权和了：

//计算4*4个整型点组成的矩形点阵的左上角点的坐标，也即(x0, y0)
int x0 = floor(x_float) - 1;
int?y0?=?floor(y_float)?-?1;


__m128?sum_m?=?_mm_setzero_ps();   //0 0 0 0
  
uchar?*src_p?=?src.ptr(y0);  //4*4矩形点阵的第一行首地址
__m128?src_m?=?_mm_set_ps(src_p[x0?+?3],?src_p[x0?+?2],?src_p[x0?+?1], src_p[x0]);?//4*4矩形点阵的第一行点像素值：A(x0+3,y0) A(x0+2,y0) A(x0+1,y0) A(x0,y0)???????????????????????
sum_m?=?_mm_add_ps(sum_m, _mm_mul_ps(src_m, coeff0));??//累加：W(3,0)*A(x0+3,y0) W(2,0)*A(x0+2,y0) W(1,0)*A(x0+1,y0) W(0,0)*A(x0,y0) 


src_p?=?src.ptr(y0?+?1);   //4*4矩形点阵的第二行首地址
src_m?=?_mm_set_ps(src_p[x0?+?3],?src_p[x0?+?2],?src_p[x0?+?1], src_p[x0]);??//4*4矩形点阵的第二行点像素值：A(x0+3,y1) A(x0+2,y1) A(x0+1,y1) A(x0,y1)
sum_m?=?_mm_add_ps(sum_m, _mm_mul_ps(src_m, coeff1)); //累加：W(3,1)*A(x0+3,y1) W(2,1)*A(x0+2,y1) W(1,1)*A(x0+1,y1) W(0,1)*A(x0,y1)


src_p?=?src.ptr(y0?+?2);   //4*4矩形点阵的第三行首地址
src_m?=?_mm_set_ps(src_p[x0?+?3],?src_p[x0?+?2],?src_p[x0?+?1], src_p[x0]);???//4*4矩形点阵的第三行点像素值：A(x0+3,y2) A(x0+2,y2) A(x0+1,y2) A(x0,y2)
sum_m?=?_mm_add_ps(sum_m, _mm_mul_ps(src_m, coeff2)); //累加：W(3,2)*A(x0+3,y2) W(2,2)*A(x0+2,y2) W(1,2)*A(x0+1,y2) W(0,2)*A(x0,y2)


src_p?=?src.ptr(y0?+?3);   //4*4矩形点阵的第四行首地址
src_m?=?_mm_set_ps(src_p[x0?+?3],?src_p[x0?+?2],?src_p[x0?+?1], src_p[x0]);???//4*4矩形点阵的第四行点像素值：A(x0+3,y3) A(x0+2,y3) A(x0+1,y3) A(x0,y3)
sum_m?=?_mm_add_ps(sum_m, _mm_mul_ps(src_m, coeff3)); //累加：W(3,3)*A(x0+3,y3) W(2,3)*A(x0+2,y3) W(1,3)*A(x0+1,y3) W(0,3)*A(x0,y3)


float?*p?=?(float?*)&sum_m;
uchar?sum?=?(uchar)(p[0]?+?p[1]?+?p[2]?+?p[3]);???//最后再把sum_m中的四个累加和加起来，即得到最终的插值结果

完整的SSE指令优化的双三次插值代码如下：

uchar cubic_inner_SSE(Mat src, float x_float, float y_float, float a)
{
  //计算权重系数
  float u = x_float - floor(x_float);
  float v = y_float - floor(y_float);
  float a_mul_4 = (a + a) + (a + a);   //4a
  float a_mul_5 = a_mul_4 + a;         //5a
  float a_mul_8 = a_mul_4 + a_mul_4;   //8a
  float a_add_3 = a + 3;
  float a_add_2 = a + 2;
  __m128 a_m = _mm_set1_ps(a);
  __m128 m_1 = _mm_set1_ps(1.0);
  __m128 a_mul_4_m = _mm_set1_ps(a_mul_4);
  __m128 a_mul_5_m = _mm_set1_ps(a_mul_5);
  __m128 a_mul_8_m = _mm_set1_ps(a_mul_8);
  __m128 a_add_3_m = _mm_set1_ps(a_add_3);
  __m128 a_add_2_m = _mm_set1_ps(a_add_2);


  __m128 C30_A30 = _mm_set_ps(2 - v, 1 + v, 2 - u, 1 + u);   //C3 C0 A3 A0
  __m128 C21_A21 = _mm_set_ps(1 - v, v, 1 - u, u);   //C2 C1 A2 A1


  __m128 tmp0 = _mm_sub_ps(_mm_mul_ps(a_m, C30_A30), a_mul_5_m);   //a*xx - a_mul_5
  tmp0 = _mm_add_ps(a_mul_8_m, _mm_mul_ps(C30_A30, tmp0));       //a_mul_8 + xx*(a*xx - a_mul_5)
  tmp0 = _mm_sub_ps(_mm_mul_ps(C30_A30, tmp0), a_mul_4_m);    //xx*(a_mul_8 + xx*(a*xx - a_mul_5)) - a_mul_4  = C3 C0 A3 A0


  __m128 tmp1 = _mm_sub_ps(_mm_mul_ps(a_add_2_m, C21_A21), a_add_3_m);   //a_add_2*xx - a_add_3
  tmp1 = _mm_mul_ps(_mm_mul_ps(C21_A21, C21_A21), tmp1);    //xx*xx*(a_add_2*xx - a_add_3)
  tmp1 = _mm_add_ps(m_1, tmp1);     //1 + xx*xx*(a_add_2*xx - a_add_3) = C2 C1 A2 A1


  __m128 A_m = _mm_unpacklo_ps(tmp0, tmp1);    //tmp1[1] tmp0[1] tmp1[0] tmp0[0] = A2 A3 A1 A0
  __m128 C_m = _mm_unpackhi_ps(tmp0, tmp1);    //tmp1[3] tmp0[3] tmp1[2] tmp0[2] = C2 C3 C1 C0
  A_m = _mm_shuffle_ps(A_m, A_m, _MM_SHUFFLE(2, 3, 1, 0));   //A3 A2 A1 A0
  C_m = _mm_shuffle_ps(C_m, C_m, _MM_SHUFFLE(2, 3, 1, 0));   //C3 C2 C1 C0


  __declspec(align(16)) float C[4];
  _mm_store_ps(C, C_m);


  __m128 m128_C = _mm_set1_ps(C[0]);
  __m128 coeff0 = _mm_mul_ps(A_m, m128_C);


  m128_C = _mm_set1_ps(C[1]);
  __m128 coeff1 = _mm_mul_ps(A_m, m128_C);


  m128_C = _mm_set1_ps(C[2]);
  __m128 coeff2 = _mm_mul_ps(A_m, m128_C);


  m128_C = _mm_set1_ps(C[3]);
??__m128?coeff3?=?_mm_mul_ps(A_m,?m128_C);
??
  ///////////////////////////////////////////////////////////////////////////////////////////////////


  int x0 = floor(x_float) - 1;
  int y0 = floor(y_float) - 1;


  __m128 sum_m = _mm_setzero_ps();
  
  uchar *src_p = src.ptr(y0);
  __m128 src_m = _mm_set_ps(src_p[x0 + 3], src_p[x0 + 2], src_p[x0 + 1], src_p[x0]);
  sum_m = _mm_add_ps(sum_m, _mm_mul_ps(src_m, coeff0));


  src_p = src.ptr(y0 + 1);
  src_m = _mm_set_ps(src_p[x0 + 3], src_p[x0 + 2], src_p[x0 + 1], src_p[x0]);
  sum_m = _mm_add_ps(sum_m, _mm_mul_ps(src_m, coeff1));


  src_p = src.ptr(y0 + 2);
  src_m = _mm_set_ps(src_p[x0 + 3], src_p[x0 + 2], src_p[x0 + 1], src_p[x0]);
  sum_m = _mm_add_ps(sum_m, _mm_mul_ps(src_m, coeff2));


  src_p = src.ptr(y0 + 3);
  src_m = _mm_set_ps(src_p[x0 + 3], src_p[x0 + 2], src_p[x0 + 1], src_p[x0]);
  sum_m = _mm_add_ps(sum_m, _mm_mul_ps(src_m, coeff3));


  float *p = (float *)&sum_m;
  uchar sum = (uchar)(p[0] + p[1] + p[2] + p[3]);


??return?sum;
}

接下来，我们分别调用以上实现的cubic_inner函数和cubic_inner_SSE函数来实现图像缩放功能，实现代码如下：

//图像缩放函数
void resize_img(Mat src, Mat &dst, float row_m, float col_m)
{
  const int row = (int)(src.rows*row_m);
  const int col = (int)(src.cols*col_m);
  const float x_a = 1.0 / col_m;
  const float y_a = 1.0 / row_m;


  Mat dst_tmp = Mat::zeros(row, col, CV_8UC1);


??for?(int?i?=?0;?i?(i);
????float?y?=?i*y_a;


??? for?(int?j?=?0;?j?

自己实现了一个微秒级计时的类，用于记录函数的运行时间：

class Timer_Us
{
  private:
    LARGE_INTEGER cpuFreq;
    LARGE_INTEGER startTime;
    LARGE_INTEGER endTime;
    
  public:
    double rumTime;
    void get_frequence(void);
    void start_timer(void);
    void stop_timer(char *str);
    Timer_Us();    //构造函数
    ~Timer_Us();   //析构函数
};


void Timer_Us::get_frequence(void)
{
  QueryPerformanceFrequency(&cpuFreq);   //获取时钟频率
}


void Timer_Us::start_timer(void)
{
  QueryPerformanceCounter(&startTime);    //开始计时
}




void Timer_Us::stop_timer(char *str)
{
  QueryPerformanceCounter(&endTime);    //结束计时
  rumTime = (((endTime.QuadPart - startTime.QuadPart) * 1000.0f) / cpuFreq.QuadPart);
  cout <

最后是测试函数，调用以上实现的图像缩放函数，对248*236的Lena图像的宽和高都放大到原来的三倍，并记录SSE指令优化插值前后的耗时。

void resize_img_test(void)
{
  Mat img = imread("lena.tif", CV_LOAD_IMAGE_GRAYSCALE);
  
  Timer_Us timer;
  float mul = 3;  //宽和高都放大三倍


  Mat img_resize;
  
  timer.start_timer();  //开始计时
  resize_img(img, img_resize2, mul, mul, 2);
  timer.stop_timer("cubic resize time:");  //结束计时，并显示运行耗时
  
  imshow("cubic img_resize", img_resize);
  waitKey();
}

运行以上代码，调用原C++实现的cubic_inner函数，耗时约35.5022 ms，如果是调用SSE指令优化的cubic_inner_SSE函数，耗时约17.3297 ms。因此SSE优化之后，耗时减少约一半，优化效果还是比较理想的。

原图

放大的图像

实际上以上实现的图像缩放函数resize_img还有很大的优化空间，比如函数里面有两层循环，外面一层是行遍历，里面一层是列遍历，在双三次插值过程中，有一些参数的计算对于同一行数据来说是一样的，因此可以把这部分计算过程从内循环提到外循环来做，如此以来，每一行只需要计算一次这些参数，可以减少不少耗时。进一步优化的resize_img函数代码如下。调用该函数对同样的Lena图像进行宽、高各三倍的放大，耗时减少为10 ms左右，优化效果还是比较显著的。

void resize_img_cubic(Mat src, Mat &dst, float row_m, float col_m)
{
  const int row = (int)(src.rows*row_m);
  const int col = (int)(src.cols*col_m);
  const float x_a = 1.0 / col_m;
  const float y_a = 1.0 / row_m;


  Mat dst_tmp = Mat::zeros(row, col, CV_8UC1);


  __declspec(align(16)) float A[4];    //内存对齐
  float C[4];
  float a = -0.15;
  //这些参数不变，直接提到循环外面计算
  float a_mul_4 = (a + a) + (a + a);   //4a
  float a_mul_5 = a_mul_4 + a;         //5a
  float a_mul_8 = a_mul_4 + a_mul_4;   //8a
  float a_add_3 = a + 3;
  float a_add_2 = a + 2;
  float xx;


  for (int i = 0; i < row; i++)
  {
    uchar *p = dst_tmp.ptr(i);
    
????//以下这些是提到外循环计算的参数
    float y = i*y_a;
    int y0 = floor(y) - 1;
????float?v?=?y?-?floor(y);
    xx = 1 + v;
    C[0] = -a_mul_4 + xx*(a_mul_8 + xx*(a*xx - a_mul_5));   //1(y0);
    uchar *src_p1 = src.ptr(y0+1);
    uchar *src_p2 = src.ptr(y0+2);
    uchar *src_p3 = src.ptr(y0+3);


    for (int j = 0; j < col; j++)
    {
      float x = j*x_a;
      float u = x - floor(x);
      
      xx = 1 + u;
      A[0] = -a_mul_4 + xx*(a_mul_8 + xx*(a*xx - a_mul_5)); //-a_mul_4 + a_mul_8*u - a_mul_5*u*u + a*u*u*u;   //1

学习代码优化有一段时间了，包括代码自身结构优化、SSE指令优化、CUDA优化等。感触最深的是，代码优化是一个精益求精的过程，一步步地优化之后，往往优化代码与原来的代码相比已经面目全非了，因此优化之后的代码可读性非常差，如果不对自己的优化思路作详细记录，过一段时间可能自己都看不懂自己的优化代码了，这是非常尴尬的，所以详细记录与注释还是非常有必要的。当然，本人的水平有限，以上代码的优化只是一个抛砖引玉的过程，也许还有更大的优化空间，如果读者有更好的优化idea，欢迎给我留言讨论。