这里的高斯模糊采用的是论文《Recursive implementation of the Gaussian filter》里描述的递归算法。
仔细观察和理解上述公式，在forward过程中,n是递增的，因此，如果在进行forward之前，把in数据先完整的赋值给w，然后式子（9a）就可以变为：
w[n] = B w[n] + (b1 w[n-1] + b2 w[n-2] + b3 w[n-3]) / b0； ---------> (1a)
在backward过程中，n是递减的，因此在backward前，把w的数据完整的拷贝到out中，则式子(9b)变为：
out[n] = B out[n] + (b1 out[n+1] + b2 out[n+2] + b3 out[n+3]) / b0 ； <--------- (1b)
从编程角度看来，backward中的拷贝是完全没有必要的，因此 (1b)可以直接写为：
w[n] = B w[n] + (b1 w[n+1] + b2 w[n+2] + b3 w[n+3]) / b0 ； <--------- (1c)
从速度上考虑，浮点除法很慢，因此，我们重新定义b1 = b1 / b0, b2 = b2 / b0, b3 = b3 / b0, 最终得到我们使用的递归公式：
w[n] = B w[n] + b1 w[n-1] + b2 w[n-2] + b3 w[n-3]； ---------> (a)
w[n] = B w[n] + b1 w[n+1] + b2 w[n+2] + b3 w[n+3] ； <--------- (b)
上述公式是一维的高斯模糊计算方法，针对二维的图像，正确的做法就是先对每个图像行进行模糊处理得到中间结果，再对中间结果的每列进行模糊操作，最终得到二维的模糊结果，当然也可以使用先列后行这样的做法。
另外注意点就是，边缘像素的处理，我们看到在公式(a)或者(b)中，w[n]的结果分别依赖于前三个或者后三个元素的值，而对于边缘位置的值，这些都是不在合理范围内的，通常的做法是镜像数据或者重复边缘数据，实践证明，由于是递归的算法，起始值的不同会将结果一直延续下去，因此，不同的方法对边缘部分的结果还是有一定的影响的，这里我们采用的简单的重复边缘像素值的方式。
由于上面公式中的系数均为浮点类型，因此，计算一般也是对浮点进行的，也就是说一般需要先把图像数据转换为浮点，然后进行高斯模糊，在将结果转换为字节类型的图像，因此，上述过程可以分别用一下几个函数完成：
CalcGaussCof // 计算高斯模糊中使用到的系数
ConvertBGR8U2BGRAF // 将字节数据转换为浮点数据
GaussBlurFromLeftToRight // 水平方向的前向传播
GaussBlurFromRightToLeft // 水平方向的反向传播
GaussBlurFromTopToBottom // 垂直方向的前向传播
GaussBlurFromBottomToTop // 垂直方向的反向传播
ConvertBGRAF2BGR8U // 将结果转换为字节数据
我们依次分析之。
CalcGaussCof，这个很简单，就按照论文中提出的计算公式根据用户输入的参数来计算，当然结合下上面我们提到的除法变乘法的优化，注意，为了后续的一些标号的问题，我讲上述公式(a)和(b)中的系数B写为b0了。
void CalcGaussCof(float Radius, float &B0, float &B1, float &B2, float &B3) { float Q, B; if (Radius >= 2.5) Q = (double)(0.98711 * Radius - 0.96330); // 对应论文公式11b else if ((Radius >= 0.5) && (Radius < 2.5)) Q = (double)(3.97156 - 4.14554 * sqrt(1 - 0.26891 * Radius)); else Q = (double)0.1147705018520355224609375; B = 1.57825 + 2.44413 * Q + 1.4281 * Q * Q + 0.422205 * Q * Q * Q; // 对应论文公式8c B1 = 2.44413 * Q + 2.85619 * Q * Q + 1.26661 * Q * Q * Q; B2 = -1.4281 * Q * Q - 1.26661 * Q * Q * Q; B3 = 0.422205 * Q * Q * Q; B0 = 1.0 - (B1 + B2 + B3) / B; B1 = B1 / B; B2 = B2 / B; B3 = B3 / B; }
由上述代码可见，B0+B1+B2+B3=1，即是归一化的系数，这也是算法可以递归进行的关键之一。
接着为了方便中间过程，我们先将字节数据转换为浮点数据，这部分代码也很简单：
void ConvertBGR8U2BGRAF(unsigned char *Src, float *Dest, int Width, int Height, int Stride) { //#pragma omp parallel for for (int Y = 0; Y < Height; Y++) { unsigned char *LinePS = Src + Y * Stride; float *LinePD = Dest + Y * Width * 3; for (int X = 0; X < Width; X++, LinePS += 3, LinePD += 3) { LinePD[0] = LinePS[0]; LinePD[1] = LinePS[1]; LinePD[2] = LinePS[2]; } } }
大家可以尝试自己把内部的X循环再展开看看效果。
水平方向的前向传播按照公式去写也是很简单的，但是直接使用公式里面有很多访问数组的过程（不一定就慢），我稍微改造下写成如下形式：
void GaussBlurFromLeftToRight(float *Data, int Width, int Height, float B0, float B1, float B2, float B3) { //#pragma omp parallel for for (int Y = 0; Y < Height; Y++) { float *LinePD = Data + Y * Width * 3; float BS1 = LinePD[0], BS2 = LinePD[0], BS3 = LinePD[0]; // 边缘处使用重复像素的方案 float GS1 = LinePD[1], GS2 = LinePD[1], GS3 = LinePD[1]; float RS1 = LinePD[2], RS2 = LinePD[2], RS3 = LinePD[2]; for (int X = 0; X < Width; X++, LinePD += 3) { LinePD[0] = LinePD[0] * B0 + BS1 * B1 + BS2 * B2 + BS3 * B3; LinePD[1] = LinePD[1] * B0 + GS1 * B1 + GS2 * B2 + GS3 * B3; // 进行顺向迭代 LinePD[2] = LinePD[2] * B0 + RS1 * B1 + RS2 * B2 + RS3 * B3; BS3 = BS2, BS2 = BS1, BS1 = LinePD[0]; GS3 = GS2, GS2 = GS1, GS1 = LinePD[1]; RS3 = RS2, RS2 = RS1, RS1 = LinePD[2]; } } }
不多描述，请大家把上述代码和公式(a)对应一下就知道了。
有了GaussBlurFromLeftToRight的参考代码，GaussBlurFromRightToLeft的代码就不会有什么大的困难了，为了避免一些懒人不看文章不思考，这里我不贴这段的代码。
那么垂直方向上简单的做只需要改变下循环的方向，以及每次指针增加量更改为Width * 3即可。
那么我们来考虑下垂直方向能不能不这样处理呢，指针每次增加Width * 3会造成严重的Cache miss，特别是对于宽度稍微大点的图像，我们仔细观察垂直方向，会发现依旧可以按照Y -- X这样的循环方式也是可以的，代码如下：
void GaussBlurFromTopToBottom(float *Data, int Width, int Height, float B0, float B1, float B2, float B3) { for (int Y = 0; Y < Height; Y++) { float *LinePD3 = Data + (Y + 0) * Width * 3; float *LinePD2 = Data + (Y + 1) * Width * 3; float *LinePD1 = Data + (Y + 2) * Width * 3; float *LinePD0 = Data + (Y + 3) * Width * 3; for (int X = 0; X < Width; X++, LinePD0 += 3, LinePD1 += 3, LinePD2 += 3, LinePD3 += 3) { LinePD0[0] = LinePD0[0] * B0 + LinePD1[0] * B1 + LinePD2[0] * B2 + LinePD3[0] * B3; LinePD0[1] = LinePD0[1] * B0 + LinePD1[1] * B1 + LinePD2[1] * B2 + LinePD3[1] * B3; LinePD0[2] = LinePD0[2] * B0 + LinePD1[2] * B1 + LinePD2[2] * B2 + LinePD3[2] * B3; } } }
就是说我们不是一下子就把列方向的前向传播进行完，而是每次只进行一个像素的传播，当一行所有像素都进行完了列方向的前向传播后，在切换到下一行，这样处理就避免了Cache miss出现的情况。
注意到列方向的边缘位置，为了方便代码的处理，我们在高度方向上上下分别扩展了3个行的像素，当进行完中间的有效行的行方向前向和反向传播后，按照前述的重复边缘像素的方式填充上下那空出的三行数据。
同样的道理，GaussBlurFromBottomToTop的代码可由读者自己补充进去。
最后的ConvertBGRAF2BGR8U也很简单，就是一个for循环：
void ConvertBGRAF2BGR8U(float *Src, unsigned char *Dest, int Width, int Height, int Stride) { //#pragma omp parallel for for (int Y = 0; Y < Height; Y++) { float *LinePS = Src + Y * Width * 3; unsigned char *LinePD = Dest + Y * Stride; for (int X = 0; X < Width; X++, LinePS += 3, LinePD += 3) { LinePD[0] = LinePS[0]; LinePD[1] = LinePS[1]; LinePD[2] = LinePS[2]; } } }
在有效的范围内，上述浮点计算的结果不会超出byte所能表达的范围，因此也不需要特别的进行Clamp操作。
最后就是一些内存分配和释放的代码了，如下所示：
void GaussBlur(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, float Radius) { float B0, B1, B2, B3; float *Buffer = (float *)malloc(Width * (Height + 6) * sizeof(float) * 3); CalcGaussCof(Radius, B0, B1, B2, B3); ConvertBGR8U2BGRAF(Src, Buffer + 3 * Width * 3, Width, Height, Stride); GaussBlurFromLeftToRight(Buffer + 3 * Width * 3, Width, Height, B0, B1, B2, B3); GaussBlurFromRightToLeft(Buffer + 3 * Width * 3, Width, Height, B0, B1, B2, B3); // 如果启用多线程，建议把这个函数写到GaussBlurFromLeftToRight的for X循环里，因为这样就可以减少线程并发时的阻力 memcpy(Buffer + 0 * Width * 3, Buffer + 3 * Width * 3, Width * 3 * sizeof(float)); memcpy(Buffer + 1 * Width * 3, Buffer + 3 * Width * 3, Width * 3 * sizeof(float)); memcpy(Buffer + 2 * Width * 3, Buffer + 3 * Width * 3, Width * 3 * sizeof(float)); GaussBlurFromTopToBottom(Buffer, Width, Height, B0, B1, B2, B3); memcpy(Buffer + (Height + 3) * Width * 3, Buffer + (Height + 2) * Width * 3, Width * 3 * sizeof(float)); memcpy(Buffer + (Height + 4) * Width * 3, Buffer + (Height + 2) * Width * 3, Width * 3 * sizeof(float)); memcpy(Buffer + (Height + 5) * Width * 3, Buffer + (Height + 2) * Width * 3, Width * 3 * sizeof(float)); GaussBlurFromBottomToTop(Buffer, Width, Height, B0, B1, B2, B3); ConvertBGRAF2BGR8U(Buffer + 3 * Width * 3, Dest, Width, Height, Stride); free(Buffer); }
正如上所述，分配了Height + 6行的内存区域，主要是为了方便垂直方向的处理，在执行GaussBlurFromTopToBottom之前按照重复边缘的原则复制3行，然后在GaussBlurFromBottomToTop之前在复制底部边缘的3行像素。
至此，一个简单而又高效的高斯模糊就基本完成了，代码数量也不多，也没有啥难度，不晓得为什么很多人搞不定。
按照我的测试，上述方式代码在一台I5-6300HQ 2.30GHZ的笔记本上针对一副3000*2000的24位图像的处理时间大约需要370ms，如果在C++的编译选项的代码生成选项里的启用增强指令集选择-->流式处理SIMD扩展2（/arch:sse2)，则编译后的程序大概需要220ms的时间。
我们尝试利用系统的一些资源进一步提高速度，首先我们想到了SSE优化，关于这方面Intel有一篇官方的文章和代码：
IIR Gaussian Blur Filter Implementation using Intel® Advanced Vector Extensions [PDF 513KB]
source code: gaussian_blur.cpp [36KB]
我只是简单的看了下这篇文章，感觉他里面用到的计算公式和Deriche滤波器的很像，和本文参考的Recursive implementation 不太一样，并且其SSE代码对能处理的图还有很多限制，对我这里的参考意义不大。
我们先看下核心的计算的SSE优化，注意到 GaussBlurFromLeftToRight 的代码中， 其核心的计算部分是几个乘法，但是他只有3个乘法计算，如果能够凑成四行，那么就可以充分利用SSE的批量计算功能了，也就是如果能增加一个通道，使得GaussBlurFromLeftToRight变为如下形式：
void GaussBlurFromLeftToRight(float *Data, int Width, int Height, float B0, float B1, float B2, float B3) { //#pragma omp parallel for for (int Y = 0; Y < Height; Y++) { float *LinePD = Data + Y * Width * 4; float BS1 = LinePD[0], BS2 = LinePD[0], BS3 = LinePD[0]; // 边缘处使用重复像素的方案 float GS1 = LinePD[1], GS2 = LinePD[1], GS3 = LinePD[1]; float RS1 = LinePD[2], RS2 = LinePD[2], RS3 = LinePD[2]; float AS1 = LinePD[3], AS2 = LinePD[3], AS3 = LinePD[3]; for (int X = 0; X < Width; X++, LinePD += 4) { LinePD[0] = LinePD[0] * B0 + BS1 * B1 + BS2 * B2 + BS3 * B3; LinePD[1] = LinePD[1] * B0 + GS1 * B1 + GS2 * B2 + GS3 * B3; // 进行顺向迭代 LinePD[2] = LinePD[2] * B0 + RS1 * B1 + RS2 * B2 + RS3 * B3; LinePD[3] = LinePD[3] * B0 + AS1 * B1 + AS2 * B2 + AS3 * B3; BS3 = BS2, BS2 = BS1, BS1 = LinePD[0]; GS3 = GS2, GS2 = GS1, GS1 = LinePD[1]; RS3 = RS2, RS2 = RS1, RS1 = LinePD[2]; AS3 = AS2, AS2 = AS1, AS1 = LinePD[3]; } } }
则很容易就把上述代码转换成SSE可以规范处理的代码了。
而对于Y方向的代码，你仔细观察会发现，无论是及通道的图，天然的就可以使用SSE进行处理，详见下面的代码。
好，我们还是一个一个的来分析：
第一个函数 CalcGaussCof 无需进行任何的优化。
第二个函数 ConvertBGR8U2BGRAF按照上述分析需要重新写，因为需要增加一个通道，新的通道的值填0或者其他值都可以，但建议填0，这对有些SSE函数很有用，我把这个函数的SSE实现共享一下：
void ConvertBGR8U2BGRAF_SSE(unsigned char *Src, float *Dest, int Width, int Height, int Stride) { const int BlockSize = 4; int Block = (Width - 2) / BlockSize; __m128i Mask = _mm_setr_epi8(0, 1, 2, -1, 3, 4, 5, -1, 6, 7, 8, -1, 9, 10, 11, -1); // Mask为-1的地方会自动设置数据为0 __m128i Zero = _mm_setzero_si128(); //#pragma omp parallel for for (int Y = 0; Y < Height; Y++) { unsigned char *LinePS = Src + Y * Stride; float *LinePD = Dest + Y * Width * 4; int X = 0; for (; X < Block * BlockSize; X += BlockSize, LinePS += BlockSize * 3, LinePD += BlockSize * 4) { __m128i SrcV = _mm_shuffle_epi8(_mm_loadu_si128((const __m128i*)LinePS), Mask); // 提取了16个字节，但是因为24位的，我们要将其变为32位的，所以只能提取出其中的前12个字节 __m128i Src16L = _mm_unpacklo_epi8(SrcV, Zero); __m128i Src16H = _mm_unpackhi_epi8(SrcV, Zero); _mm_store_ps(LinePD + 0, _mm_cvtepi32_ps(_mm_unpacklo_epi16(Src16L, Zero))); // 分配内存时已经是16字节对齐了，然后每行又是4的倍数的浮点数，因此，每行都是16字节对齐的 _mm_store_ps(LinePD + 4, _mm_cvtepi32_ps(_mm_unpackhi_epi16(Src16L, Zero))); // _mm_stream_ps是否快点？ _mm_store_ps(LinePD + 8, _mm_cvtepi32_ps(_mm_unpacklo_epi16(Src16H, Zero))); _mm_store_ps(LinePD + 12, _mm_cvtepi32_ps(_mm_unpackhi_epi16(Src16H, Zero))); } for (; X < Width; X++, LinePS += 3, LinePD += 4) { LinePD[0] = LinePS[0]; LinePD[1] = LinePS[1]; LinePD[2] = LinePS[2]; LinePD[3] = 0; // Alpha最好设置为0，虽然在下面的CofB0等SSE常量中通过设置ALPHA对应的系数为0，可以使得计算结果也为0，但是不是最合理的 } } }
稍作解释，_mm_loadu_si128一次性加载16个字节的数据到SSE寄存器中，对于24位图像，16个字节里包含了5 + 1 / 3个像素的信息，而我们的目标是把这些数据转换为4通道的信息，因此，我们只能一次性的提取处16/4=4个像素的值，这借助于_mm_shuffle_epi8函数和合适的Mask来实现，_mm_unpacklo_epi8/_mm_unpackhi_epi8分别提取了SrcV的高8位和低8位的8个字节数据并将它们转换为8个16进制数保存到Src16L和Src16H中， 而_mm_unpacklo_epi16/_mm_unpackhi_epi16则进一步把16位数据扩展到32位整形数据，最后通过_mm_cvtepi32_ps函数把整形数据转换为浮点型。
可能有人注意到了 int Block = (Width - 2) / BlockSize; 这一行代码，为什么要-2操作呢，这也是我在多次测试发现程序总是出现内存错误时才意识到的，因为_mm_loadu_si128一次性加载了5 + 1 / 3个像素的信息，当在处理最后一行像素时（其他行不会），如果Block 取Width/BlockSize, 则很有可能访问了超出像素范围内的内存，而-2不是-1就是因为那个额外的1/3像素起的作用。
接着看水平方向的前向传播的SSE方案：
void GaussBlurFromLeftToRight_SSE(float *Data, int Width, int Height, float B0, float B1, float B2, float B3) { const __m128 CofB0 = _mm_set_ps(0, B0, B0, B0); const __m128 CofB1 = _mm_set_ps(0, B1, B1, B1); const __m128 CofB2 = _mm_set_ps(0, B2, B2, B2); const __m128 CofB3 = _mm_set_ps(0, B3, B3, B3); //#pragma omp parallel for for (int Y = 0; Y < Height; Y++) { float *LinePD = Data + Y * Width * 4; __m128 V1 = _mm_set_ps(LinePD[3], LinePD[2], LinePD[1], LinePD[0]); __m128 V2 = V1, V3 = V1; for (int X = 0; X < Width; X++, LinePD += 4) // 还有一种写法不需要这种V3/V2/V1递归赋值，一次性计算3个值，详见 D:程序设计正在研究CoreShockFilterConvert里的高斯模糊，但速度上没啥区别 { __m128 V0 = _mm_load_ps(LinePD); __m128 V01 = _mm_add_ps(_mm_mul_ps(CofB0, V0), _mm_mul_ps(CofB1, V1)); __m128 V23 = _mm_add_ps(_mm_mul_ps(CofB2, V2), _mm_mul_ps(CofB3, V3)); __m128 V = _mm_add_ps(V01, V23); V3 = V2; V2 = V1; V1 = V; _mm_store_ps(LinePD, V); } } }
和上面的4通道的GaussBlurFromLeftToRight_SSE比较，你会发现基本上语法上没有任何的区别，实在是太简单了，注意我没有用_mm_storeu_ps，而是直接使用_mm_store_ps，这是因为，第一，分配Data内存时，我采用了_mm_malloc分配了16字节对齐的内存，而Data每行元素的个数又都是4的倍数，因此，每行起始位置处的内存也是16字节对齐的，所以直接_mm_store_ps完全是可以的。
同理，在垂直方向的前向传播的SSE优化代码就更直接了：
void GaussBlurFromTopToBottom_SSE(float *Data, int Width, int Height, float B0, float B1, float B2, float B3) { const __m128 CofB0 = _mm_set_ps(0, B0, B0, B0); const __m128 CofB1 = _mm_set_ps(0, B1, B1, B1); const __m128 CofB2 = _mm_set_ps(0, B2, B2, B2); const __m128 CofB3 = _mm_set_ps(0, B3, B3, B3); for (int Y = 0; Y < Height; Y++) { float *LinePS3 = Data + (Y + 0) * Width * 4; float *LinePS2 = Data + (Y + 1) * Width * 4; float *LinePS1 = Data + (Y + 2) * Width * 4; float *LinePS0 = Data + (Y + 3) * Width * 4; for (int X = 0; X < Width * 4; X += 4) { __m128 V3 = _mm_load_ps(LinePS3 + X); __m128 V2 = _mm_load_ps(LinePS2 + X); __m128 V1 = _mm_load_ps(LinePS1 + X); __m128 V0 = _mm_load_ps(LinePS0 + X); __m128 V01 = _mm_add_ps(_mm_mul_ps(CofB0, V0), _mm_mul_ps(CofB1, V1)); __m128 V23 = _mm_add_ps(_mm_mul_ps(CofB2, V2), _mm_mul_ps(CofB3, V3)); _mm_store_ps(LinePS0 + X, _mm_add_ps(V01, V23)); } } }
对上面的代码不想做任何解释了。
最有难度的应该算是ConvertBGRAF2BGR8U的SSE版本了，由于某些原因，我不太愿意分享这个函数的代码，有兴趣的朋友可以参考opencv的有关实现。
最后的SSE版本高斯模糊的主要代码如下：
void GaussBlur_SSE(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, float Radius) { float B0, B1, B2, B3; float *Buffer = (float *)_mm_malloc(Width * (Height + 6) * sizeof(float) * 4, 16); CalcGaussCof(Radius, B0, B1, B2, B3); ConvertBGR8U2BGRAF_SSE(Src, Buffer + 3 * Width * 4, Width, Height, Stride); GaussBlurFromLeftToRight_SSE(Buffer + 3 * Width * 4, Width, Height, B0, B1, B2, B3); // 在SSE版本中，这两个函数占用的时间比下面两个要多,不过C语言版本也是一样的 GaussBlurFromRightToLeft_SSE(Buffer + 3 * Width * 4, Width, Height, B0, B1, B2, B3); // 如果启用多线程，建议把这个函数写到GaussBlurFromLeftToRight的for X循环里，因为这样就可以减少线程并发时的阻力 memcpy(Buffer + 0 * Width * 4, Buffer + 3 * Width * 4, Width * 4 * sizeof(float)); memcpy(Buffer + 1 * Width * 4, Buffer + 3 * Width * 4, Width * 4 * sizeof(float)); memcpy(Buffer + 2 * Width * 4, Buffer + 3 * Width * 4, Width * 4 * sizeof(float)); GaussBlurFromTopToBottom_SSE(Buffer, Width, Height, B0, B1, B2, B3); memcpy(Buffer + (Height + 3) * Width * 4, Buffer + (Height + 2) * Width * 4, Width * 4 * sizeof(float)); memcpy(Buffer + (Height + 4) * Width * 4, Buffer + (Height + 2) * Width * 4, Width * 4 * sizeof(float)); memcpy(Buffer + (Height + 5) * Width * 4, Buffer + (Height + 2) * Width * 4, Width * 4 * sizeof(float)); GaussBlurFromBottomToTop_SSE(Buffer, Width, Height, B0, B1, B2, B3); ConvertBGRAF2BGR8U_SSE(Buffer + 3 * Width * 4, Dest, Width, Height, Stride); _mm_free(Buffer); }
对于同样的3000*2000的彩色图像，SSE版本的代码耗时只有145ms的耗时，相对于普通的C代码有约2.5倍左右的提速，这也情有可原，因为我们在执行SSE的代码时时多处理了一个通道的计算量的，但是同编译器自身的SSE优化220ms，只有1.5倍的提速了，这说明编译器的SSE优化能力还是相当强的。
进一步的优化就是我上面的注释掉的opemp相关代码，在ConvertBGR8U2BGRAF / GaussBlurFromLeftToRight / GaussBlurFromRightToLeft / ConvertBGRAF2BGR8U 4个函数中，直接加上简单的#pragma omp parallel for就可以了，至于GaussBlurFromTopToBottom_SSE、 GaussBlurFromBottomToTop_SSE则由于上下行之间数据的相关性，是无法实现并行计算的，但是测试表示他们的耗时要比水平方向的少很多。
比如我们指定openmp使用2个线程，在上述机器上（四核的)，纯C版本能优化到252ms，而纯SSE的只能提高到100ms左右，编译器自身的SSE优化速度大约是150ms，基本上还是保持同一个级别的比例。
对于灰度图像，很可惜，上述的水平方向上的优化方式就无论为力了，一种方式就是把4行灰度像素混洗成一行，但是这个操作不太方便用SSE实现，另外一种就是把水平方向的数据先转置，然后利用垂直方向的SSE优化代码处理，完成在转置回去，最后对转置的数据再次进行垂直方向SSE优化，当然转置的过程是可以借助于SSE的代码实现的，但是需要额外的一份内存，速度上可能和普通的C相比就不会有那么多的提升了，这个待有时间了再去测试。
前前后后写这个大概也花了半个月的时间，写完了整个流程，在倒过来看，真的是非常的简单，有的时候就是这样。
本文并没有提供完整的可以直接执行的全部代码，需者自勤，提供一段C#的调用工程供有兴趣的朋友测试和比对（未使用OPENMP版本）。
http://files.cnblogs.com/files/Imageshop/GaussBLur_Sample.rar
后记：后来测试发现，当半径参数较大时，无论是C版本还是SSE版本都会出现一些不规则的瑕疵，感觉像是溢出了，后来发现主要是当半径变大后，B0参数变得很小，以至于用float类型的数据来处理递归已经无法保证足够的精度了，解决的办法是使用double类型，这对于C语言版本来说是秒秒的事情，而对于SSE版本则是较大的灾难，double时换成AVX可能改动量不大，但是AVX的普及度以及.....，不过，一般情况下，半径不大于75时结果都还是正确的，这对于大部分的应用来说是足够的，半径75时，整个图像已经差不多没有任何的细节了，再大，区别也不大了。
解决上述问题一个可行的方案就是使用Deriche滤波器，我用该滤波器的float版本对大半径是不会出现问题的，并且也有相关SSE参考代码。
后续文章：高斯模糊算法的全面优化过程分享（二）。