这儿的高斯含糊选用的是论文《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:\程序设计\正在研讨\Core\ShockFilter\Convert里的高斯含糊，但速度上没啥差异
{
__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://share.eepw.com.cn/share/download/id/383730
跋文：后来测验发现，当半径参数较大时，不管是C版别仍是SSE版别都会呈现一些不规则的瑕疵，感觉像是溢出了，后来发现首要是当半径变大后，B0参数变得很小，以至于用float类型的数据来处理递归现已无法确保满足的精度了，处理的方法是运用double类型，这关于C言语版别来说是秒秒的工作，而关于SSE版别则是较大的灾祸，double时换成AVX或许改动量不大，可是AVX的遍及度以及…..，不过，一般状况下，半径不大于75时成果都仍是正确的，这关于大部分的运用来说是满足的，半径75时，整个图画现已差不多没有任何的细节了，再大，差异也不大了。
处理上述问题一个可行的计划便是运用Deriche滤波器，我用该滤波器的float版别对大半径是不会呈现问题的，而且也有相关SSE参阅代码。
