SSE优化系列十一：SSE优化三次卷积插值算法

在阅读了https://www.cnblogs.com/Imageshop/p/9069650.html 博主的三次卷积插值的进一步SSE优化文章后&＃xff0c;对这个算法产生了很大的兴趣&＃xff0c;然而复现代码的过程是艰难的&＃xff0c;好在作者了提供了部分代码&＃xff0c;经过1周的努力&＃xff0c;我实现了完整的支持1通道&＃xff0c;3通道&＃xff0c;4通道的三次卷积插值算法的SSE优化代码。现在准备分享一下这个算法的实现过程。

之前我在我的另外一个工程: https://github.com/BBuf/Image-processing-algorithm/tree/master/Image%20Interpolation 实现过3次卷积插值&＃xff0c;百度百科提供的原理如下&＃xff1a;https://baike.baidu.com/item/%E5%8F%8C%E4%B8%89%E6%AC%A1%E6%8F%92%E5%80%BC/11055947?fr&＃61;aladdin 。这个过程可以用另外一套理论来解释&＃xff0c;也就是&＃xff1a;

图片截图自ImageShop的文章。先贴出使用三次卷积插值的简单代码&＃xff1a;

void Bicubic_Original(unsigned char *Src, int Width, int Height, int Stride, unsigned char *Pixel, float X, float Y)
{int Channel &＃61; Stride / Width;int PosX &＃61; floor(X), PosY &＃61; floor(Y);float PartXX &＃61; X - PosX, PartYY &＃61; Y - PosY;unsigned char *Pixel00 &＃61; GetCheckedPixel(Src, Width, Height, Stride, Channel, PosX - 1, PosY - 1);unsigned char *Pixel01 &＃61; GetCheckedPixel(Src, Width, Height, Stride, Channel, PosX &＃43; 0, PosY - 1);unsigned char *Pixel02 &＃61; GetCheckedPixel(Src, Width, Height, Stride, Channel, PosX &＃43; 1, PosY - 1);unsigned char *Pixel03 &＃61; GetCheckedPixel(Src, Width, Height, Stride, Channel, PosX &＃43; 2, PosY - 1);unsigned char *Pixel10 &＃61; GetCheckedPixel(Src, Width, Height, Stride, Channel, PosX - 1, PosY &＃43; 0);unsigned char *Pixel11 &＃61; GetCheckedPixel(Src, Width, Height, Stride, Channel, PosX &＃43; 0, PosY &＃43; 0);unsigned char *Pixel12 &＃61; GetCheckedPixel(Src, Width, Height, Stride, Channel, PosX &＃43; 1, PosY &＃43; 0);unsigned char *Pixel13 &＃61; GetCheckedPixel(Src, Width, Height, Stride, Channel, PosX &＃43; 2, PosY &＃43; 0);unsigned char *Pixel20 &＃61; GetCheckedPixel(Src, Width, Height, Stride, Channel, PosX - 1, PosY &＃43; 1);unsigned char *Pixel21 &＃61; GetCheckedPixel(Src, Width, Height, Stride, Channel, PosX &＃43; 0, PosY &＃43; 1);unsigned char *Pixel22 &＃61; GetCheckedPixel(Src, Width, Height, Stride, Channel, PosX &＃43; 1, PosY &＃43; 1);unsigned char *Pixel23 &＃61; GetCheckedPixel(Src, Width, Height, Stride, Channel, PosX &＃43; 2, PosY &＃43; 1);unsigned char *Pixel30 &＃61; GetCheckedPixel(Src, Width, Height, Stride, Channel, PosX - 1, PosY &＃43; 2);unsigned char *Pixel31 &＃61; GetCheckedPixel(Src, Width, Height, Stride, Channel, PosX &＃43; 0, PosY &＃43; 2);unsigned char *Pixel32 &＃61; GetCheckedPixel(Src, Width, Height, Stride, Channel, PosX &＃43; 1, PosY &＃43; 2);unsigned char *Pixel33 &＃61; GetCheckedPixel(Src, Width, Height, Stride, Channel, PosX &＃43; 2, PosY &＃43; 2);float U0 &＃61; SinXDivX(1 &＃43; PartXX), U1 &＃61; SinXDivX(PartXX);float U2 &＃61; SinXDivX(1 - PartXX), U3 &＃61; SinXDivX(2 - PartXX);float V0 &＃61; SinXDivX(1 &＃43; PartYY), V1 &＃61; SinXDivX(PartYY);float V2 &＃61; SinXDivX(1 - PartYY), V3 &＃61; SinXDivX(2 - PartYY);for (int I &＃61; 0; I }


从这里可以看到&＃xff0c;我们引入上面那2条曲线正是为了计算每个像素点在三次卷积插值时的系数。而提出一个拟合得表达式来近似Sin曲线的原因是什么呢&＃xff1f;实际上当我们将X取值为0.3&＃xff0c;然后按照这2个公式计算的结果如下&＃xff1a;
SinXDivX_Standard(1 &＃43; X) &＃43; SinXDivX_Standard(X) &＃43; SinXDivX_Standard(1 - X) &＃43; SinXDivX_Standard(2 - X) &＃61; 0.8767
SinXDivX(1 &＃43; X) &＃43; SinXDivX(X) &＃43; SinXDivX(1 - X) &＃43; SinXDivX(2 - X) &＃61;1&＃xff0c; 可以看到如果我们使用拟合得表达式&＃xff0c;权重系数就不需要再进行归一化处理了。这两个函数的代码如下&＃xff1a;

// 该函数计算插值曲线sin(x * PI) / (x * PI)的值,下面是它的近似拟合表达式
float SinXDivX(float X) {const float a &＃61; -1; //a还可以取 a&＃61;-2,-1,-0.75,-0.5等等&＃xff0c;起到调节锐化或模糊程度的作用X &＃61; abs(X);float X2 &＃61; X * X, X3 &＃61; X2 * X;if (X <&＃61; 1)return (a &＃43; 2) * X3 - (a &＃43; 3) * X2 &＃43; 1;else if (X <&＃61; 2)return a * X3 - (5 * a) * X2 &＃43; (8 * a) * X - (4 * a);elsereturn 0;
}// 精确计算插值曲线sin(x * PI) / (x * PI)
float SinXDivX_Standard(float X) {if (abs(X) <0.000001f)return 1;elsereturn sin(X * 3.1415926f) / (X * 3.1415926f);
}

三次卷积插值的原始实现

// 原始的插值算法
void IM_Resize_Cubic_Origin(unsigned char *Src, unsigned char *Dest, int SrcW, int SrcH, int StrideS, int DstW, int DstH, int StrideD) {int Channel &＃61; StrideS / SrcW;if ((SrcW &＃61;&＃61; DstW) && (SrcH &＃61;&＃61; DstH)) {memcpy(Dest, Src, SrcW * SrcH * Channel * sizeof(unsigned char));return;}printf("%d\n", Channel);for (int Y &＃61; 0; Y }

优化1&＃xff1a;边界特殊处理&＃43;查表法

对原始实现一个显然的优化点在于&＃xff0c;对于每个像素都要计算领域每个像素的系数&＃xff0c;显然这个可以用定点化查表来实现。同时&＃xff0c;如果我们把边界特殊处理&＃xff0c;我们就可以省去很多GetCheckedPixel函数的调用&＃xff0c;这个加速也是非常明显的。最后在使用了边界特殊处理&＃43;定点化查表后&＃xff0c;速度比原始实现大概提升了2倍。

// C语言实现的查表&＃43;插值算法
void IM_Resize_Cubic_Table(unsigned char *Src, unsigned char *Dest, int SrcW, int SrcH, int StrideS, int DstW, int DstH, int StrideD) {int Channel &＃61; StrideS / SrcW;if ((SrcW &＃61;&＃61; DstW) && (SrcH &＃61;&＃61; DstH)) {memcpy(Dest, Src, SrcW * SrcH * Channel * sizeof(unsigned char));return;}short *SinXDivX_Table &＃61; (short *)malloc(513 * sizeof(short));for (int I &＃61; 0; I <513; I&＃43;&＃43;)SinXDivX_Table[I] &＃61; int(0.5 &＃43; 256 * SinXDivX(I / 256.0f)); // 建立查找表&＃xff0c;定点化int AddX &＃61; (SrcW <<16) / DstW, AddY &＃61; (SrcH <<16) / DstH;int ErrorX &＃61; -(1 <<15) &＃43; (AddX >> 1), ErrorY &＃61; -(1 <<15) &＃43; (AddY >> 1);int StartX &＃61; ((1 <<16) - ErrorX) / AddX &＃43; 1; // 计算出需要特殊处理的边界int StartY &＃61; ((1 <<16) - ErrorY) / AddY &＃43; 1; // y0&＃43;y*yr>&＃61;1; y0&＃61;ErrorY &＃61;> y>&＃61;(1-ErrorY)/yrint EndX &＃61; (((SrcW - 3) <<16) - ErrorX) / AddX &＃43; 1;int EndY &＃61; (((SrcH - 3) <<16) - ErrorY) / AddY &＃43; 1; // y0&＃43;y*yr<&＃61;(height-3) &＃61;> y<&＃61;(height-3-ErrorY)/yrif (StartY >&＃61; DstH) StartY &＃61; DstH;if (StartX >&＃61; DstW) StartX &＃61; DstW;if (EndX }

优化3&＃xff1a;SSE优化

如果是一个通道和三个通道的图像&＃xff0c;是很容易把上面的过程翻译为SSE代码的。比较复杂的是对于24位图像的处理&＃xff0c;我这里偷了个懒&＃xff0c;将RGB图像直接转为了RGBA图像&＃xff0c;在RGBA图像操作之后将结果图像又转为了RGB图像&＃xff0c;这个转化过程也使用了SSE优化。在SSE优化系列2-高斯滤波https://blog.csdn.net/just_sort/article/details/95212099中&＃xff0c;留下了一个坑&＃xff0c;就是BGR和BGRA相互转化的SSE优化&＃xff0c;我这里填了这个坑点。实现的代码如下&＃xff1a;

void ConvertBGR8U2BGRAF_SSE(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride) {const int BlockSize &＃61; 4;int Block &＃61; (Width - 2) / BlockSize;__m128i Mask &＃61; _mm_setr_epi8(0, 1, 2, -1, 3, 4, 5, -1, 6, 7, 8, -1, 9, 10, 11, -1);__m128i Mask2 &＃61; _mm_setr_epi8(0, 2, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1);__m128i Zero &＃61; _mm_setzero_si128();for (int Y &＃61; 0; Y }void ConvertBGRAF2BGR8U_SSE(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride) {const int BlockSize &＃61; 4;int Block &＃61; (Width - 2) / BlockSize;//__m128i Mask &＃61; _mm_setr_epi8(0, 1, 2, 4, 5, 6, 8, 9, 10, 12, 13, 14, 3, 7, 11, 15);__m128i MaskB &＃61; _mm_setr_epi8(0, 4, 8, 12, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1);__m128i MaskG &＃61; _mm_setr_epi8(1, 5, 9, 13, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1);__m128i MaskR &＃61; _mm_setr_epi8(2, 6, 10, 14, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1);__m128i Zero &＃61; _mm_setzero_si128();for (int Y &＃61; 0; Y }


都是一些常规操作&＃xff0c;可以看着数据模拟一下寄存器变量变化过程。这里再提供一个输出寄存器变量值的函数&＃xff1a;

void debug(__m128i var) {uint8_t *val &＃61; (uint8_t*)&var;//can also use uint32_t instead of 16_t printf("Numerical: %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i %i\n",val[0], val[1], val[2], val[3], val[4], val[5],val[6], val[7], val[8], val[9], val[10], val[11], val[12], val[13],val[14], val[15]);
}


最后对单通道和四通道的图像进行三次卷积插值的SSE优化代码为&＃xff1a;

// 使用SSE优化立方插值算法
// 最大支持图像大小为: 32767*32767
void IM_Resize_SSE(unsigned char *Src, unsigned char *Dest, int SrcW, int SrcH, int StrideS, int DstW, int DstH, int StrideD) {int Channel &＃61; StrideS / SrcW;if ((SrcW &＃61;&＃61; DstW) && (SrcH &＃61;&＃61; DstH)) {memcpy(Dest, Src, SrcW * SrcH * Channel * sizeof(unsigned char));return;}short *SinXDivX_Table &＃61; (short *)malloc(513 * sizeof(short));short *Table &＃61; (short *)malloc(DstW * 4 * sizeof(short));for (int I &＃61; 0; I <513; I&＃43;&＃43;)SinXDivX_Table[I] &＃61; int(0.5 &＃43; 256 * SinXDivX(I / 256.0f)); // 建立查找表&＃xff0c;定点化int AddX &＃61; (SrcW <<16) / DstW, AddY &＃61; (SrcH <<16) / DstH;int ErrorX &＃61; -(1 <<15) &＃43; (AddX >> 1), ErrorY &＃61; -(1 <<15) &＃43; (AddY >> 1);int StartX &＃61; ((1 <<16) - ErrorX) / AddX &＃43; 1; // 计算出需要特殊处理的边界int StartY &＃61; ((1 <<16) - ErrorY) / AddY &＃43; 1; // y0&＃43;y*yr>&＃61;1; y0&＃61;ErrorY &＃61;> y>&＃61;(1-ErrorY)/yrint EndX &＃61; (((SrcW - 3) <<16) - ErrorX) / AddX &＃43; 1;int EndY &＃61; (((SrcH - 3) <<16) - ErrorY) / AddY &＃43; 1; // y0&＃43;y*yr<&＃61;(height-3) &＃61;> y<&＃61;(height-3-ErrorY)/yrif (StartY >&＃61; DstH) StartY &＃61; DstH;if (StartX >&＃61; DstW) StartX &＃61; DstW;if (EndX > 8);Table[X * 4 &＃43; 0] &＃61; SinXDivX_Table[256 &＃43; U]; //建立一个新表便于SSE操作Table[X * 4 &＃43; 1] &＃61; SinXDivX_Table[U];Table[X * 4 &＃43; 2] &＃61; SinXDivX_Table[256 - U];Table[X * 4 &＃43; 3] &＃61; SinXDivX_Table[512 - U];}int SrcY &＃61; ErrorY;for (int Y &＃61; 0; Y > 8);unsigned char *LineY &＃61; Src &＃43; ((SrcY >> 16) - 1) * StrideS;__m128i PartY &＃61; _mm_setr_epi32(SinXDivX_Table[256 &＃43; V], SinXDivX_Table[V], SinXDivX_Table[256 - V], SinXDivX_Table[512 - V]);for (int X &＃61; StartX; X > 16) - 1) * Channel;unsigned char *Pixel1 &＃61; Pixel0 &＃43; StrideS;unsigned char *Pixel2 &＃61; Pixel1 &＃43; StrideS;unsigned char *Pixel3 &＃61; Pixel2 &＃43; StrideS;if (Channel &＃61;&＃61; 1) {__m128i P01 &＃61; _mm_cvtepu8_epi16(_mm_unpacklo_epi32(_mm_cvtsi32_si128(*((int *)Pixel0)), _mm_cvtsi32_si128(*((int *)Pixel1)))); // P00 P01 P02 P03 P10 P11 P12 P13__m128i P23 &＃61; _mm_cvtepu8_epi16(_mm_unpacklo_epi32(_mm_cvtsi32_si128(*((int *)Pixel2)), _mm_cvtsi32_si128(*((int *)Pixel3)))); // P20 P21 P22 P23 P30 P31 P32 P33__m128i Sum01 &＃61; _mm_madd_epi16(P01, PartX); // P00 * U0 &＃43; P01 * U1 P02 * U2 &＃43; P03 * U3 P10 * U0 &＃43; P11 * U1 P12 * U2 &＃43; P13 * U3__m128i Sum23 &＃61; _mm_madd_epi16(P23, PartX); // P20 * U0 &＃43; P21 * U1 P22 * U2 &＃43; P23 * U3 P30 * U0 &＃43; P31 * U1 P32 * U2 &＃43; P33 * U3__m128i Sum &＃61; _mm_hadd_epi32(Sum01, Sum23); // P00 * U0 &＃43; P01 * U1 &＃43; P02 * U2 &＃43; P03 * U3 P10 * U0 &＃43; P11 * U1 &＃43; P12 * U2 &＃43; P13 * U3 P20 * U0 &＃43; P21 * U1 &＃43; P22 * U2 &＃43; P23 * U3 P30 * U0 &＃43; P31 * U1 &＃43; P32 * U2 &＃43; P33 * U3LinePD[0] &＃61; ClampToByte(_mm_hsum_epi32(_mm_mullo_epi32(Sum, PartY)) >> 16);}else if (Channel &＃61;&＃61; 4) {__m128i P0 &＃61; _mm_loadu_si128((__m128i *)Pixel0), P1 &＃61; _mm_loadu_si128((__m128i *)Pixel1);__m128i P2 &＃61; _mm_loadu_si128((__m128i *)Pixel2), P3 &＃61; _mm_loadu_si128((__m128i *)Pixel3);P0 &＃61; _mm_shuffle_epi8(P0, _mm_setr_epi8(0, 4, 8, 12, 1, 5, 9, 13, 2, 6, 10, 14, 3, 7, 11, 15)); // B0 G0 R0 A0P1 &＃61; _mm_shuffle_epi8(P1, _mm_setr_epi8(0, 4, 8, 12, 1, 5, 9, 13, 2, 6, 10, 14, 3, 7, 11, 15)); // B1 G1 R1 A1P2 &＃61; _mm_shuffle_epi8(P2, _mm_setr_epi8(0, 4, 8, 12, 1, 5, 9, 13, 2, 6, 10, 14, 3, 7, 11, 15)); // B2 G2 R2 A2P3 &＃61; _mm_shuffle_epi8(P3, _mm_setr_epi8(0, 4, 8, 12, 1, 5, 9, 13, 2, 6, 10, 14, 3, 7, 11, 15)); // B3 G3 R3 A3__m128i BG01 &＃61; _mm_unpacklo_epi32(P0, P1); // B0 B1 G0 G1__m128i RA01 &＃61; _mm_unpackhi_epi32(P0, P1); // R0 R1 A0 A1__m128i BG23 &＃61; _mm_unpacklo_epi32(P2, P3); // B2 B3 G2 G3__m128i RA23 &＃61; _mm_unpackhi_epi32(P2, P3); // R2 R3 A2 A3__m128i B01 &＃61; _mm_unpacklo_epi8(BG01, _mm_setzero_si128());__m128i B23 &＃61; _mm_unpacklo_epi8(BG23, _mm_setzero_si128());__m128i SumB &＃61; _mm_hadd_epi32(_mm_madd_epi16(B01, PartX), _mm_madd_epi16(B23, PartX));__m128i G01 &＃61; _mm_unpackhi_epi8(BG01, _mm_setzero_si128());__m128i G23 &＃61; _mm_unpackhi_epi8(BG23, _mm_setzero_si128());__m128i SumG &＃61; _mm_hadd_epi32(_mm_madd_epi16(G01, PartX), _mm_madd_epi16(G23, PartX));__m128i R01 &＃61; _mm_unpacklo_epi8(RA01, _mm_setzero_si128());__m128i R23 &＃61; _mm_unpacklo_epi8(RA23, _mm_setzero_si128());__m128i SumR &＃61; _mm_hadd_epi32(_mm_madd_epi16(R01, PartX), _mm_madd_epi16(R23, PartX));__m128i A01 &＃61; _mm_unpackhi_epi8(RA01, _mm_setzero_si128());__m128i A23 &＃61; _mm_unpackhi_epi8(RA23, _mm_setzero_si128());__m128i SumA &＃61; _mm_hadd_epi32(_mm_madd_epi16(A01, PartX), _mm_madd_epi16(A23, PartX));__m128i Result &＃61; _mm_setr_epi32(_mm_hsum_epi32(_mm_mullo_epi32(SumB, PartY)), _mm_hsum_epi32(_mm_mullo_epi32(SumG, PartY)), _mm_hsum_epi32(_mm_mullo_epi32(SumR, PartY)), _mm_hsum_epi32(_mm_mullo_epi32(SumA, PartY)));Result &＃61; _mm_srai_epi32(Result, 16);// *((int *)LinePD) &＃61; _mm_cvtsi128_si32(_mm_packus_epi16(_mm_packus_epi32(Result, Result), Result));_mm_stream_si32((int *)LinePD, _mm_cvtsi128_si32(_mm_packus_epi16(_mm_packus_epi32(Result, Result), Result)));//LinePD[0] &＃61; IM_ClampToByte(_mm_hsum_epi32(_mm_mullo_epi32(SumB, PartY)) >> 16); // 确实有部分存在超出unsigned char范围的&＃xff0c;因为定点化的缘故//LinePD[1] &＃61; IM_ClampToByte(_mm_hsum_epi32(_mm_mullo_epi32(SumG, PartY)) >> 16);//LinePD[2] &＃61; IM_ClampToByte(_mm_hsum_epi32(_mm_mullo_epi32(SumR, PartY)) >> 16);//LinePD[3] &＃61; IM_ClampToByte(_mm_hsum_epi32(_mm_mullo_epi32(SumA, PartY)) >> 16);}}for (int X &＃61; EndX; X }


这个优化过程技巧性比较强&＃xff0c;这里需要仔细说明一下。我们先来看Channel为1的情况&＃xff0c;观察这句代码&＃xff1a;Pixel00[I] * U0 &＃43; Pixel01[I] * U1 &＃43; Pixel02[I] * U2 &＃43; Pixel03[I] * U3&＃xff0c;我们发现Pixel00,Pixel01,Pixel02,Pixel03在内存中是连续的&＃xff0c;且这个变量都是在0-256的值域中&＃xff0c;显然我们可以使用SSE寄存器读取这几个变量。同时SSE中有一个_mm_madd_epi16&＃xff0c;实现的功能是&＃xff1a;

__m128i _mm_madd_epi16 (__m128i a, __m128i b); Return ValueAdds the signed 32-bit integer results pairwise and packs the 4 signed 32-bit integer results.r0 :&＃61; (a0 * b0) &＃43; (a1 * b1)r1 :&＃61; (a2 * b2) &＃43; (a3 * b3)r2 :&＃61; (a4 * b4) &＃43; (a5 * b5)r3 :&＃61; (a6 * b6) &＃43; (a7 * b7)


考虑我们函数中&＃xff0c;行1到行4每行的代码都只有4次乘法和3次加法&＃xff0c;不能直接使用&＃xff0c;但是我们可以考虑把两行整合在一起&＃xff0c;一次性计算&＃xff0c;这样就需要调用2次_mm_madd_epi16 &＃xff0c;然后2次的结果在调用_mm_hadd_epi32这个水平方向的累加函数就能得到新的结果&＃xff0c; _mm_hsum_epi32的实现如下:

// 4个有符号的32位的数据相加的和
inline int _mm_hsum_epi32(__m128i V) { //V3 V2 V1 V0__m128i T &＃61; _mm_add_epi32(V, _mm_srli_si128(V, 8)); //V3&＃43;V1 V2&＃43;V0 V1 V0T &＃61; _mm_add_epi32(T, _mm_srli_si128(T, 4)); //V3&＃43;V1&＃43;V2&＃43;V0 V2&＃43;V0&＃43;V1 V1&＃43;V0 V0return _mm_cvtsi128_si32(T); //提取低位
}


实现过程中还有很多以前重复介绍过的sse指令&＃xff0c;可以在MSDN或我的资源中&＃xff1a; https://github.com/BBuf/Image-processing-algorithm-Speed/tree/master/resources 查看。算法的具体过程请看代码&＃xff0c; 一些重要地方我都注释好了。完整代码请在github查看。

可以看到我们的SSE优化版本代码还是比不过OpenCV3.1.0自带的函数&＃xff0c;猜测OpenCV内部对BGR图像应该是直接在3个通道上操作&＃xff0c;且优化做得更多更好&＃xff0c;所以比我SSE优化的速度快了近3倍。

原图&＃xff1a;

结果图&＃xff0c;扩大了1.5倍&＃xff1a;

https://www.cnblogs.com/Imageshop/p/9069650.html