【滤波·5】超高速指数模糊算法

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

1. 算法核心及原理

德里奇过滤器 (tsdconseil.fr)

两个链接的核心代码就是

    *zR += (alpha * ((R << zprec) - *zR)) >> aprec;

事实上可以看出代码实现的就是 Deriche filter里提到的公式：  a first order IIR, which equation is: y[n] = a*x[n] + (1-a)*y[n-1]  。

一阶的alpha参数代码里使用的是：

int alpha = (int)((1 << aprec) * (1.0f - std::exp(-2.3f / (radius + 1.f))));

而恰好系数也是一直以来我存疑的地方，不知道该怎么设，自己使用的是

总结：  也就是指数模糊实际上就是derich版的迭代高斯滤波，只不过实现的是特定alpha系数的一阶公式，同时使用了定点化，利用zprec，aprec参数控制精度。

2.实现

/*
其中Pixel就是我们要处理的像素，zR,zG,zB,zA是外部传入的一个累加值，alpha、aprec、zprec是由模糊半径radius生成的一些固定的系数。

　　类似于高斯模糊或者StackBlur，算法也是属于行列分离的，行方向上，先用上述方式从左向右计算，然后在从右向左，接着进行列方向处理，先从上向下，然后在从下向上。当然，行列的计算也可以反过来。需要注意的是，每一步都是用之前处理过的数据进行的。
*/
static inline void _blurinner(guchar* pixel, gint* zR,gint* zG,gint* zB, gint* zA, gint alpha,gint aprec,gint zprec)
{
	gint R, G, B, A;
	R = *pixel;
	G = *(pixel + 1);
	B = *(pixel + 2);
	A = *(pixel + 3);

	*zR += (alpha * ((R << zprec) - *zR)) >> aprec;
	*zG += (alpha * ((G << zprec) - *zG)) >> aprec;
	*zB += (alpha * ((B << zprec) - *zB)) >> aprec;
	*zA += (alpha * ((A << zprec) - *zA)) >> aprec;

	*pixel       = *zR >> zprec;
	*(pixel + 1) = *zG >> zprec;
	*(pixel + 2) = *zB >> zprec;
	*(pixel + 3) = *zA >> zprec;
}
  
  
static inline void _blurrow(guchar* pixels,
                             gint    width,
                             gint    /* height */,  // TODO: This seems very strange. Why is height not used as it is in _blurcol() ?
                             gint    channels,
                             gint    line,
                             gint    alpha,
                             gint    aprec,
                             gint    zprec)
{
	gint    zR;
	gint    zG;
	gint    zB;
	gint    zA;
	gint    index;
	guchar* scanline;

	scanline = &(pixels[line * width * channels]);

	zR = *scanline << zprec;
	zG = *(scanline + 1) << zprec;
	zB = *(scanline + 2) << zprec;
	zA = *(scanline + 3) << zprec;

	for (index = 0; index < width; index ++)
	_blurinner(&scanline[index * channels], &zR, &zG, &zB, &zA, alpha, aprec,
	zprec);

	for (index = width - 2; index >= 0; index--)
	_blurinner(&scanline[index * channels], &zR, &zG, &zB, &zA, alpha, aprec,
	zprec);
}

static inline void _blurcol(guchar* pixels,
                               gint    width,
                               gint    height,
                               gint    channels,
                               gint    x,
                               gint    alpha,
                               gint    aprec,
                               gint    zprec)
{
	gint zR;
	gint zG;
	gint zB;
	gint zA;
	gint index;
	guchar* ptr;

	ptr = pixels;

	ptr += x * channels;

	zR = *((guchar*) ptr    ) << zprec;
	zG = *((guchar*) ptr + 1) << zprec;
	zB = *((guchar*) ptr + 2) << zprec;
	zA = *((guchar*) ptr + 3) << zprec;

	for (index = width; index < (height - 1) * width; index += width)
	  _blurinner((guchar*) &ptr[index * channels], &zR, &zG, &zB, &zA, alpha,
	  aprec, zprec);

	for (index = (height - 2) * width; index >= 0; index -= width)
	  _blurinner((guchar*) &ptr[index * channels], &zR, &zG, &zB, &zA, alpha,
	  aprec, zprec);
}

// pixels   image-data
// width    image-width
// height   image-height
// channels image-channels
// in-place blur of image 'img' with kernel of approximate radius 'radius'
// blurs with two sided exponential impulse response
// aprec = precision of alpha parameter in fixed-point format 0.aprec
// zprec = precision of state parameters zR,zG,zB and zA in fp format 8.zprecb  
void _expblur(guchar* pixels,
                 gint    width,
                 gint    height,
                 gint    channels,
                 gint    radius,
                 gint    aprec,
                 gint    zprec)
{
	gint alpha;
	gint row = 0;
	gint col = 0;

	if (radius < 1)
	  return;

	// calculate the alpha such that 90% of 
	// the kernel is within the radius.
	// (Kernel extends to infinity)
	alpha = (gint) ((1 << aprec) * (1.0f - expf(-2.3f / (radius + 1.f))));
	//常用的aprec取值为16，Zprec 取值为7。

	for (; row < height; row++)
	  _blurrow(pixels, width, height, channels, row, alpha, aprec, zprec);

	for (; col < width; col++)
	  _blurcol(pixels, width, height, channels, col, alpha, aprec, zprec);

	return;
}

作为一个典型的应用，或者说尽量减少参数，常用的aprec取值为16，Zprec 取值为7。

     回顾下代码，整体过程中除了alpha参数的计算涉及到了浮点，其他部分都是整形的乘法和移位操作，因此可以想象，速度应该不慢，而且非常适合于手机端处理。同时注意到_blurrow和_blurcol函数循环明显相互之间是独立的，可以利用多线程并行处理，但是这个代码主要是专注于算法的表达，并没有过多的考虑更好的效率。

     另外一点，很明显，算法的耗时是和Radius参数没有任何关系的，也就是说这也是个O(1)算法。

3.灰度图处理

3.1水平方向

for (int Y = 0; Y < Height; Y++)
{
    byte *LinePS = Src + Y * Stride;
    byte *LinePD = Dest + Y * Stride;
    int Sum = LinePS[0] << Zprec;
    for (int X = 0; X < Width; X++)　　　　　　//　　从左往右
    {
        Sum += (Alpha * ((LinePS[X] << Zprec) - Sum)) >> Aprec;
        LinePD[X] = Sum >> Zprec;
    }
    for (int X = Width - 1; X >= 0; X--)　　　//　　从右到左
    {
        Sum += (Alpha * ((LinePD[X] << Zprec) - Sum)) >> Aprec;
        LinePD[X] = Sum >> Zprec;
    }
}

3.2垂直方向

在 高斯模糊算法的全面优化过程分享（一） 中我们探讨过垂直方向处理算法一般不宜直接写，而应该用一个临时的行缓存进行处理，这样列方向的灰度图的处理代码类似下面的：

int *Buffer = (int *)malloc(Width * sizeof(int));
for (int X = 0; X < Width; X++)        Buffer[X] = Src[X] << Zprec;
for (int Y = 0; Y < Height; Y++)
{
    byte *LinePS = Src + Y * Stride;
    byte *LinePD = Dest + Y * Stride;
    for (int X = 0; X < Width; X++)　　　　　　　　//　　从上到下
    {
        Buffer[X] += (Alpha * ((LinePS[X] << Zprec) - Buffer[X])) >> Aprec;
        LinePD[X] = Buffer[X] >> Zprec;
    }
}
for (int Y = Height - 1; Y >= 0; Y--)　　　　　　//　　从下到上
{
    byte *LinePD = Dest + Y * Stride;
    for (int X = 0; X < Width; X++)
    {
        Buffer[X] += (Alpha * ((LinePD[X] << Zprec) - Buffer[X])) >> Aprec;
        LinePD[X] = Buffer[X] >> Zprec;
    }
}
free(Buffer);

       修改为上述后，测试一个3000*2000的8位灰度图，耗时大约52ms（未使用多线程的），和普通的C语言实现的Boxblur时间差不多。

       除了线程外，这个时间是否还有改进的空间呢，我们先来看看列方向的优化。

　　 在列方向的  for (int X = 0; X < Width; X++) 循环内，我们注意到针对Buffer的每个元素的处理都是独立和相同的，很明显这样的过程是很容易使用SIMD指令优化的，但是循环体内部有一些是unsigned char类型的数据，为使用SIMD指令，需要转换为int类型较为方便，而最后保存时又需要重新处理为unsigned char类型的，这种来回转换的耗时和其他计算的提速能否来带来效益呢，我们进行了代码的编写，比如：

for (int X = 0; X < Width; X++)　　　　　　　　//　　从上到下
{
    Buffer[X] += (Alpha * ((LinePS[X] << Zprec) - Buffer[X])) >> Aprec;
    LinePD[X] = Buffer[X] >> Zprec;
}

这段代码可以用如下的SIMD指令代替：

int X = 0;
for (X = 0; X < Width - 8; X += 8)            
{
    //    将8个字节数存入到2个XMM寄存器中
    //    方案1：使用SSE4新增的_mm_cvtepu8_epi32的函数，优点是两行是独立的
    __m128i Dst1 = _mm_cvtepu8_epi32(_mm_cvtsi32_si128((*(int *)(LinePD + X + 0))));    //    _mm_cvtsi32_si128把参数中的32位整形数据移到XMM寄存器的最低32位，其他为清0。
    __m128i Dst2 = _mm_cvtepu8_epi32(_mm_cvtsi32_si128((*(int *)(LinePD + X + 4))));    //    _mm_cvtepu8_epi32将低32位的整形数的4个字节直接扩展到XMM的4个32位中去。
    __m128i Buf1 = _mm_loadu_si128((__m128i *)(Buffer + X + 0));
    __m128i Buf2 = _mm_loadu_si128((__m128i *)(Buffer + X + 4));
    Buf1 = _mm_add_epi32(_mm_srai_epi32(_mm_mullo_epi32(_mm_sub_epi32(_mm_slli_epi32(Dst1, Zprec), Buf1), Alpha128), Aprec), Buf1);
    Buf2 = _mm_add_epi32(_mm_srai_epi32(_mm_mullo_epi32(_mm_sub_epi32(_mm_slli_epi32(Dst2, Zprec), Buf2), Alpha128), Aprec), Buf2);
    _mm_storeu_si128((__m128i *)(Buffer + X + 0), Buf1);
    _mm_storeu_si128((__m128i *)(Buffer + X + 4), Buf2);
    _mm_storel_epi64((__m128i *)(LinePD + X), _mm_packus_epi16(_mm_packs_epi32(_mm_srai_epi32(Buf1, Zprec), _mm_srai_epi32(Buf2, Zprec)), Zero));
}
for (; X < Width; X++)        
{
    Buffer[X] += (Alpha * ((LinePD[X] << Zprec) - Buffer[X])) >> Aprec;
    LinePD[X] = Buffer[X] >> Zprec;
}

原来的三四行代码一下子变成了几十行的代码，会不会变慢呢，其实不用担心，SIMD真的很强大，测试的结果是3000*2000的图耗时降低到42ms左右，而且垂直方向的耗时占比有原先的60%降低到了35%左右，现在的核心就是水平方向的耗时了。

     当图像不是灰度模式时，对于垂直方向的处理和灰度不会有区别，这是因为，只需要增加循环的长度就可以了。

参考

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