SSE3

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

在上篇⽂章中，我们讲解了常见的最邻近插值算法、双线性插值算法和双三次插值算法的原理与实现，三种插值算法中双三次插值算法的插值效果最

好，但其也是三种算法中计算复杂度最⾼、耗时最长的算法。本⽂在给出双三次插值C++代码的基础上，着重讲解如何使⽤SSE指令来优化该算法，

并使⽤双三次插值来实现图像的缩放，⽐较SSE指令优化前后的耗时。

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

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

floatcubic_w_f(floatx,floata)

{

if(x<=1)

{

return1-(a+3)*x*x+(a+2)*x*x*x;

}

elif(x<2)

{

return-4*a+8*a*x-5*a*x*x+a*x*x*x;

}

return0.0;

}

voidcal_cubic_coeff(floatx,floaty,float*coeff)

{

floatu=x-floor(x);

floatv=y-floor(y);

u+=1;

v+=1;

floata=-0.15;

floatA[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(ints=0;s<4;s++)

{

floatC=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;

}

}

ucharcubic_inner(Matsrc,floatx_float,floaty_float,floata)

{

floatcoeff[16];

cal_cubic_coeff(x_float,y_float,coeff);//计算权重系数

floatsum=0.0;

intx0=floor(x_float)-1;

inty0=floor(y_float)-1;

for(inti=0;i<4;i++)

{

for(intj=0;j<4;j++)

{

sum+=coeff[i*4+j]*(y0+i)[x0+j];

}

}

return((uchar)sum);

}

指令优化算法

⾸先，我们来看⼀下浮点型坐标点周围的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指令并⾏计

算代码如下：

floatu=x_float-floor(x_float);//计算x坐标额⼩数部分

floatv=y_float-floor(y_float);//计算y坐标额⼩数部分

floata_mul_4=(a+a)+(a+a);//提前计算权重公式中的4a

floata_mul_5=a_mul_4+a;//提前计算权重公式中的5a

floata_mul_8=a_mul_4+a_mul_4;//提前计算权重公式中的8a

floata_add_3=a+3;//提前计算权重公式中的a+3

floata_add_2=a+2;//提前计算权重公式中的a+2

__m128a_m=_mm_t1_ps(a);//aaaa

__m128m_1=_mm_t1_ps(1.0);//1.01.01.01.0

__m128a_mul_4_m=_mm_t1_ps(a_mul_4);//4a4a4a4a

__m128a_mul_5_m=_mm_t1_ps(a_mul_5);//5a5a5a5a

__m128a_mul_8_m=_mm_t1_ps(a_mul_8);//8a8a8a8a

__m128a_add_3_m=_mm_t1_ps(a_add_3);//a+3a+3a+3a+3

__m128a_add_2_m=_mm_t1_ps(a_add_2);//a+2a+2a+2a+2

__m128C30_A30=_mm_t_ps(2-v,1+v,2-u,1+u);//dy3dy0dx3dx0

__m128C21_A21=_mm_t_ps(1-v,v,1-u,u);//dy2dy1dx2dx1

__m128tmp0=_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)

__m128tmp1=_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进⾏重新打包与排列：

__m128A_m=_mm_unpacklo_ps(tmp0,tmp1);//交替打包tmp0与tmp1的低位数据：tmp1[1]tmp0[1]tmp1[0]tmp0[0]=w(dx2)w(dx3)w(dx1)w(dx0)

__m128C_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))floatC[4];

_mm_store_ps(C,C_m);//w(dy3)w(dy2)w(dy1)w(dy0)

__m128m128_C=_mm_t1_ps(C[0]);//w(dy0)w(dy0)w(dy0)w(dy0)

__m128coeff0=_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_t1_ps(C[1]);//w(dy1)w(dy1)w(dy1)w(dy1)

__m128coeff1=_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_t1_ps(C[2]);//w(dy2)w(dy2)w(dy2)w(dy2)

__m128coeff2=_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_t1_ps(C[3]);//w(dy3)w(dy3)w(dy3)w(dy3)

__m128coeff3=_mm_mul_ps(A_m,m128_C);//w(dx3)*w(dy3)w(dx2)*w(dy3)w(dx1)*w(dy3)w(dx0)*w(d英语四级考试查询 y3)

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

//计算4*4个整型点组成的矩形点阵的左上⾓点的坐标，也即(x0,y0)

intx0=floor(x_float)-1;

inty0=floor(y_float)-1;

__m128sum_m=_mm_tzero_ps();//0000

uchar*src_p=(y0);//4*4矩形点阵的第⼀⾏⾸地址

__m128src_m=恢的组词 _mm_t_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=(y0+1);//4*4矩形点阵的第⼆⾏⾸地址

src_m=_mm_t_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=(y0+2);//4*4矩形点阵的第三⾏⾸地址

src_m=_mm_t_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=(y0+3);//4*4矩形点阵的第四⾏⾸地址

src_m=_mm_t_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;

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

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

ucharcubic_inner_SSE(Matsrc,floatx_float,floaty_float,floata)

{

//计算权重系数

floatu=x_float-floor(x_float);

floatv=y_float-floor(y_float);

floata_mul_4=(a+a)+(a+a);//4a

floata_mul_5=a_mul_4+a;//5a

floata_mul_8=a_mul_4+a_mul_4;//8a

floata_add_3=a+3;

floata_add_2=a+2;

__m128a_m=_mm_t1_ps(a);

__m128m_1=_mm_t1_ps(1.0);

__m128a_mul_4_m=_mm_t1_ps(a_mul_4);

__m128a_mul_5_m=_mm_t1_ps(a_mul_5);

__m128a_mul_8_m=_mm_t1_ps(a_mul_8);

__m128a_add_3_m=_mm_t1_ps(a_add_3);

__m128a_add_2_m=_mm_t1_ps(a_add_2);

__m128C30_A30=_mm_t_ps(2-v,1+v,2-u,1+u);//C3C0A3A0

__m128C21_A21=_mm_t_ps(1-v,v,1-u,u);//C2C1A2A1

__m128tmp0=_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=C3C0A3A0

__m128tmp1=_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)=C2C1A2A1

__m128A_m=_mm_unpacklo_ps(tmp0,tmp1);//tmp1[1]tmp0[1]tmp1[0]tmp0[0]=A2A3A1A0

__m128C_m=_mm_unpackhi_ps(tmp0,tmp1);//tmp1[3]tmp0[3]tmp1[2]tmp0[2]=C2C3C1C0

A_m=_mm_shuffle_ps(A_m,A_m,_MM_SHUFFLE(2,3,1,0));//A3A2A1A0

C驼色配什么颜色好看 _m=_mm_shuffle_ps(C_m,C_m,_MM_SHUFFLE(2,3,1,0));//C3C2C1C0

__declspec(align(16))floatC[4];

_mm_store_ps(C,C_m);

__m128m128_C=_mm_t1_ps(C[0]);

__m128coeff0=_mm_mul_ps(A_m,m128_C);

m128_C=_mm_t1_ps(C[1]);

__m128coeff1=_mm_mul_ps(A_m,m128_C);

m128_C=_mm_t1_ps(C[2]);

__m128coeff2=_mm_mul_ps(A_m,m128_C);

m128_C=_mm_t1_ps(C[3]);

__m128coeff3=_mm_mul_ps(A_m,m128_C);

///

intx0=floor(x_float)-1;

inty0=floor(y_float)-1;

__m128sum_m=_mm_tzero_ps();

uchar*src_p=(y0);

__m128src_m=_mm_t_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=(y0+1);

src_m=_mm_t_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=(y0+2);

src_m=_mm_t_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=(y0+3);

src_m=_mm_t_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;

ucharsum=(uchar)(p[0]+p[1]+p[2]+p[3]);

returnsum;

}

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

//图像缩放函数

voidresize_img(Matsrc,Mat&dst,floatrow_m,floatcol_m)

{

constintrow=(int)(*row_m);

constintcol=(int)(*col_m);

constfloatx_a=1.0/col_m;

constfloaty_a=1.0/row_m;

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

for(inti=0;i

{

uchar*p=dst_(i);

floaty=i*y_a;

for(intj=0;j

{

floatx=j*x_a;

//p[j]=cubic_inner(src,x,y,-0.5);//原函数

p[j]=cubic_inner_SSE(src,x,y,-0.5);//SSE优化函数

}

}

dst_(dst);

}

⾃⼰实现了⼀个微秒级计时的类，⽤于记录函数的运⾏时间：

classTimer_Us

{

private:

LARGE_INTEGERcpuFreq;

LARGE_INTEGERstartTime;

LARGE_INTEGERendTime;

public:

doublerumTime;

voidget_frequence(void);

voidstart_timer(void);

voidstop_timer(char*str);

Timer_Us();//构造函数

~Timer_Us();//析构函数

};

voidTimer_Us::get_frequence(void)

{

QueryPerformanceFrequency(&cpuFreq);//获取时钟频率

}

voidTimer_Us::start_timer(void)

{

QueryPerformanceCounter(&startTime);//开始计时

}

voidTimer_Us::stop_timer(char*str)

{

QueryPerformanceCounter(&endTime);//结束计时

rumTime=(((rt)*1000.0f)/rt);

cout<

}

Timer_Us::Timer_Us()//构造函数

{

QueryPerformanceFrequency(&cpuFreq);

}

Timer_Us::~Timer_Us()//析构函数

{

}

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

voidresize_im祝寿贺词70岁老人 g_test(void)

{

Matimg=imread("",CV_LOAD_IMAGE_GRAYSCALE);

Timer_Ustimer;

floatmul=3;//宽和⾼都放⼤三倍

Matimg_resize;

_timer();//开始计时

resize_img(img,img_resize2,mul,mul,2);

_timer("cubicresizetime:");//结束计时，并显⽰运⾏耗时

imshow("cubicimg_resize",img_resize);

waitKey();

}

运⾏以上代码，调⽤原C++实现的cubic_inn历史上的海兰珠 er函数干炸里脊 ，耗时约35.5022ms，如果是调⽤SSE指令优化的cubic_inner_SSE函数，耗时约17.3297ms。

因此SSE优化之后，耗时减少约⼀半，优化效果还是⽐较理想的。

原图

放⼤的图像

实际上以上实现的图像缩放函数resize_img还有很⼤的优化空间，⽐如函数⾥⾯有两层循环，外⾯⼀层是⾏遍历，⾥⾯⼀层是列遍历，在双三次插值

过程中，有⼀些参数的计算对于同⼀⾏数据来说是⼀样的，因此可以把这部分计算过程从内循环提到外循环来做，如此以来，每⼀⾏只需要计算⼀次

这些参数，可以减少不少耗时。进⼀步优化的resize_img函数代码如下。调⽤该函数对同样的Lena图像进⾏宽、⾼各三倍的放⼤，耗时减少为10ms

左右，优化效果还是⽐较显著的。

voidresize_img_cubic(Matsrc,Mat&dst,floatrow_m,floatcol_m)

{

constintrow=(int)(*row_m);

constintcol=(int)(*col_m);

constfloatx_a=1.0/col_m;

constfloaty_a=1.0/row_m;

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

__declspec(align(16))floatA[4];//内存对齐

floatC[4];

floata=-0.15;

//这些参数不变，直接提到循环外⾯计算

floata_mul_4=(a+a)+(a+a);//4a

floata_mul_5=a_mul_4+a;//5a

floata_mul_8=a_mul_4+a_mul_4;//8a

floata_add_3=a+3;

floata_add_2=a+2;

floatxx;

for(inti=0;i

{

uchar*p=dst_(i);

//以下这些是提到外循环计算的参数

floaty=i*y_a;

inty0=floor(y)-1;

floatv=y-floor(y);

xx=1+v;

C[0]=-a_mul_4+xx*(a_mul_8+xx*(a*xx-a_mul_5));//1

xx=v;//0

C[1]迫不及待造句 =1+xx*xx*(a_add_2*xx-a_add_3);

xx=1-v;//0

C[2]=1+xx*xx*(a_add_2*xx-a_add_3);

xx=2-v;//1

C[3]=-a_mul_4+xx*(a_mul_8+xx*(a*xx-a_mul_5));

__m128m128_C0=_mm_t1_ps(C[0]);

__m128m128_C1=_mm_t1_ps(C[1]);

__m128m128_C2=_mm_t1_ps(C[2]);

__m128m128_C3=_mm_t1_ps(C[3]);

uchar*src_p0=(y0);

uchar*src_p1=(y0+1);

uchar*src_p2=(y0+2);

uchar*src_p3=(y0+3);

for(intj=0;j

{

floatx=j*x_a;

floatu=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

xx=u;//0

A[1]=1+xx*xx*(a_add_2*xx-a_add_3);//1-a_add_3*xx*xx+a_add_2*xx*xx*xx;

xx=1-u;//0

A[2]=1+xx*xx*(a_add_2*xx-a_add_3);//1-a_add_3*xx*xx+a_add_2*xx*xx*xx;

xx=2-u;//1

A[3]=-a_mul_4+xx*(a_mul_8+xx*(a*xx-a_mul_5));//-a_mul_4+a_mul_8*xx-a_mul_5*xx*xx+a*xx*xx*xx;

__m128m128_A=_mm_load_ps(A);

__m128coeff0=_mm_mul_ps(m128_A,m128_C0);

__m128coeff1=_mm_mul_ps(m128_A,m128_C1);

__m128coeff2=_mm_mul_ps(m128_A,m128_C2);

__m128coeff3=_mm_mul_ps(m128_A,m128_C3);

intx0=floor(x)-1;

__m128src_m=_mm_t_ps(src_p0[x0+3],src_p0[x0+2],src_p0[x0+1],src_p0[x0]);

__m128sum_m=_mm_mul_ps(src_m,coeff0);

src_m=_mm_t_ps(src_p1[x0+3],src_p1[x0+2],src_p1[x0+1],src_p1[x0]);

sum_m=_mm_add_ps(sum_m,_mm_mul_ps(src_m,coeff1));

src_m=_mm_t_ps(src_p2[x0+3],src_p2[x0+2],src_p2[x0+1],src_p2[x0]);

sum_m=_mm_add_ps(sum_m,_mm_mul_ps(src_m,coeff2));

src_m=_mm_t_ps(src_p3[x0+3],src_p3[x0+2],src_p3[x0+1],src_p3[x0]);

sum_m=_mm_add_ps(sum_m,_mm_mul_ps(src_m,coeff3));

float*p1=(float*)&sum_m;

p[j]=(uchar)(p1[0]+p1[1]+p1[2]+p1[3]);

}

}

dst_(dst);

}

学习代码优化有⼀段时间了，包括代码⾃⾝结构优化、SSE指令优化、CUDA优化等。感触最深的是，代码优化是⼀个精益求精的过程，⼀步步地优

化之后，往往优化代码与原来的代码相⽐已经⾯⽬全⾮了，因此优化之后的代码可读性⾮常差，如果不对⾃⼰的优化思路作详细记录，过⼀段时间可

能⾃⼰都看不懂⾃⼰的优化代码了，这是⾮常尴尬的，所以详细记录与注释还是⾮常有必要的。当然，本⼈的⽔平有限，以上代码的优化只是⼀个抛

砖引⽟的过程，也许还有更⼤的优化空间，如果读者有更好的优化idea，欢迎给我留⾔讨论。

微信公众号如下，欢迎扫码关注，欢迎私信技术交流：