新闻中心

EEPW首页 > 嵌入式系统 > 牛人业话 > 高斯模糊算法的全面优化过程分享(一)

高斯模糊算法的全面优化过程分享(一)

作者:时间:2017-03-09来源:网络

  这里的模糊采用的是论文《Recursive implementation of the Gaussian filter》里描述的递归算法。

本文引用地址:http://www.eepw.com.cn/article/201703/345016.htm

  仔细观察和理解上述公式,在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参考代码。



关键词: 高斯 模糊算法

评论

技术专区

关闭