SSE图像算法优化系列三：超高速导向滤波实现过程纪要（欢迎挑战）
自从何凯明提出导向滤波后，因为其算法的简单性和有效性，该算法得到了广泛的应用，以至于新版的matlab都将其作为标准自带的函数之一了，利用他可以解决的所有的保边滤波器的能解决的问题，比如细节增强、HDR压缩、细节羽化、去雾、风格化，而且由于其保边特性，如果很多传统函数中使用高斯滤波或者均值滤波的地方用他代替，能很好解决一些强边缘的过渡不自然问题，比如retinex、Highlight/shadow等应用中，因此，快速的实现该算法具有很强的适用意义。
本文简要的记录了本人在优化导向滤波实现的过程中所适用的优化方式和一些细节，以免时间久了后自己都不记得了，但是请不要向我直接索取源代码。
自认为目前我优化的速度在CPU版本中很难有人能超越了（仅仅使用CPU、不用多线程，下采样率0.2），如果谁有更快的算法，在第三方公证的情况下，我愿意提供1000元奖励^_^。
何凯明在导向滤波一文的相关资料中提供了其matlab代码，或者用下面的流程也可以清晰的表达：
我们看到了上面的6次取mean计算的过程，也就是浮点数的boxfilter，这个东西已经是老掉牙的一个算法了，我在几年前研究过opencv内部的这个算法，并且提出了一种比opencv实现更快的方法，详见解析opencv中Box Filter的实现并提出进一步加速的方案（源码共享） 一文。不过那里的处理时针对字节数据的，其内部用到了一些整形数据的SSE优化，如果原始数据是浮点数，那反而就更为简易了，因为SSE指令生来就是为浮点数服务的。
但是即使是这样，由于6次计算以及中间的其他一些浮点运算，依然给整个算法带来了很大的运算开销和内存开销，在很多场合还是无法满足需求的，比如实时去雾等场景。在早期我的快速去雾实现中，都是先利用下采样图的导向滤波结果，然后再双线性插值放大得到大图的透射率图，虽然在视觉效果上能解决去雾算法的速度问题，但是如果是其他场景的导向滤波需求，还是会看到很多瑕疵的。
何凯明在2015又发表了一篇《Fast Guided Filter》的文章，阐述了一种很实用的更快速的导向滤波流程：
我刚刚提的在去雾中我实用的小Trick实际上就是第六步及第七步不同，我的方式可表达如下：
6： q = meana. * I + meanb
7: q = fupsample(q, s)
很明显，由于I的参与计算，何的做法能更大程度上保持结果和原汁原味的类似，而我的方式则会产生较大的块状相似，所以人家大神就是大神。
在何的论文中已经说明下采样比例 s 取4时，计算的结果和准确结果也还是非常靠近的，我在我的实现里s 取到了5。
这样改动后，所有的boxfilter均是对下取样后的数据进行处理，当s=4时，计算量减少到原有的1/16，而s=5，则减少到了原有的1/25,当然这个时候多了2个下取样和2个上取样的算法，下取样由于是缩小，计算量很小，无需关注，而上采样，计算量和原图大小有关，根据我的测评，这个上采样的耗时可能占整个过程的一般时间左右的，是非常值得注意优化的。
首先，第一，步骤6中的两个采样过程不要分开写，直接写到同一个for循环内部，这样可以节省很多坐标的计算过程，第二，这里一般的上采样通常采用双线性插值就OK了，网络上有很多关于双线性插值的SSE优化的代码，但是那些基本都是针对32位的图像做的优化，搬到24位和8位中是不适用的，而我们会在50%以上的概率中遇到24位图像，所以说啊，网络上的东西虽多，但精华太少。
我采用的一个优化方式时，先进行水平方向的上采样到一个缓冲区中（Width * SmallH）,然后在从这个缓冲区中沿着高度方向缓冲到（Width * Height)，如下图所示：
-----------------> 
-----------------> 
由于这个上采样是针对浮点型的数据，所以中间的精度损失问题可以不用考虑，而如果是图像的字节数据，则要慎重了。
由上面的第一个图到第二个图的大概代码如下：
for (int Y = 0; Y < SmallH; Y++)
{
float PosX = 0;
float AddX = (SmallW - 1.0f) / Width; // 主要是为了减少下面插值向右增1的指针超过范围，但这样做其实是和精确的算法有一点点差异的
float *LinePDA = TempA + Y * Width * Channel; // TempA和TempB为临时分配的大小为（SmallH * Width * Channel * sizeof(float）大小的内存
float *LinePDB = TempB + Y * Width * Channel;
float *LinePA = MeanA + Y * SmallW * Channel;
float *LinePB = MeanB + Y * SmallW * Channel;
if (Channel == 1)
{
for (int X = 0; X < Width; X++)
{
int XX = (int)PosX;
float PartXX = PosX - XX;
float InvertXX = 1 - PartXX;
float *PtLeftA = LinePA + XX;
float *PtLeftB = LinePB + XX;
LinePDA[X] = PtLeftA[0] * InvertXX + PtLeftA[1] * PartXX;
LinePDB[X] = PtLeftB[0] * InvertXX + PtLeftB[1] * PartXX;
PosX += AddX;
}
}
// ...................
}
这段代码用SSE去优化的伤害的脑细胞有点多，而且由于其计算量不是太大，意义可能有限。
而由第二个图到第三个图的过程大概可有用下面的代码表述：
for (int Y = 0; Y < Height; Y++)
{
float PosY = Y * (SmallH - 1.0f) / Height;
int YY = (int)PosY;
float PartYY = PosY - YY;
float InvertYY = 1 - PartYY;
byte *LinePS = Guide + Y * Stride;
byte *LinePD = Dest + Y * Stride;
float *PtTopA = TempA + YY * Width * Channel;
float *PtBottomA = PtTopA + Width * Channel;
float *PtTopB = TempB + YY * Width * Channel;
float *PtBottomB = PtTopB + Width * Channel;
for (int X = 0;; X < Width * Channel; X++)
{
float ValueA = PtTopA[X] * InvertYY + PtBottomA[X] * PartYY;
float ValueB = PtTopB[X] * InvertYY + PtBottomB[X] * PartYY;
LinePD[X] = IM_ClampFHtoByte(ValueA * LinePS[X] + ValueB * 255);
}
}
注意最后的IM_ClampFHtoByte函数是将括号内的值限制在0和255之间的。
有很多朋友可能不知道，如果把上面的IM_ClampFHtoByte这个函数去掉，直接使用括号内的代码，VS的编译器可以很好的对上面代码进行向量化编译（VS编译只要你没有把代码生成--》启用增强指令集设置成无增强指令/arch:IA32，哪怕设置为未设置，都会把浮点的代码编译为SIMD相关指令的），而如果我们对不同的Channel，比如3通道4通道在循环里展开后，很不幸，按照我们的展开循环的理论，速度应该加快，但事实却相反了。所以我们需要充分掌握编译器的向量化特性，就能写成更高效的代码。
由于在计算过程中确实存在一些结果超出了0和255的范围，因此如果把IM_ClampFHtoByte函数去除，对有些图像会出现噪点，因此，我们不能完全依赖编译器的向量化优化了，这就必须自己写SIMD指令，由于SIMD自带了饱和处理的相关函数，而上述内部的X 的for循环是很容易用SSE处理的，唯一需要注意的就是需要把LinePS对应的字节数据转换为浮点数据，这里我简单的提示可以用如下指令将8个字节数据转换为8个浮点数：
__m128i SrcI = _mm_unpacklo_epi8(_mm_loadl_epi64((__m128i const *)(LinePS + X)), Zero); // Load the lower 64 bits of the value pointed to by p into the lower 64 bits of the result, zeroing the upper 64 bits of the result.
__m128 SrcFL = _mm_cvtepi32_ps(_mm_unpacklo_epi16(SrcI, Zero)); // 转换为浮点
__m128 SrcFH = _mm_cvtepi32_ps(_mm_unpackhi_epi16(SrcI, Zero));
里面的浮点计算的过程的SSE代码就和普通的函数调用没什么却别，最后的写到LinePD这个字节数据的过程可以用_mm_storel_epi64以及有关移位搞定。
这里这样做的另外一个好处是在Y循环中计算是独立的，因此都可以使用OPENMP加速。
使用SSE优化能将上述过程提速2倍以上。
另外一个问题，在上面的流程2的第一步中，对boxfilter的半径r也是进行了同比例的缩小的，注意到boxfilter的半径通常情况下我们都是用的整数，如果缩小后的r'也进行取整的话，举例来说，对于s =4的情况下，半径为8、9、10、11这四种情况最终得到的导向滤波结果就完全一样了，似乎这不符合我们对算法严谨性的要求，所以我们要支持一种浮点半径的boxfilter。
普通意义的boxfilter肯定是无法支持浮点半径的（这不同于高斯模糊），一种变通的办法就是取浮点半径前后的两个整形半径值做模糊，然后再线性插值，举个例子，如果下取样后的半径为4.4，则分别计算R1 = boxfilter(4)以及R2 = boxfilter(5)，最后合成得到结果R:
R = R1 * (1 - 0.4) + R2 * 0.4;
如此处理后，在大部分情况下（除了下取样后的半径为整数，比如原有半径为12，s=4，这是r'=3），计算量又会稍微增加一点，需要计算小图的12次boxfilter了，不过何必纠结这个了呢。
关于上述浮点版本的Boxfilter，其实还有一种更好的实现方式。我在13行代码实现最快速最高效的积分图像算法中也提供了一段实现方框模糊的代码，当然那个代码还不是最优的，因为其中的pixlecount需要每个像素都重新计算，其实当半径较小时中间部分的像素的pixlecount为固定值，因此可以把边缘部分的像素特殊处理，对于本例，是需要进行的浮点版本的算法，那对于中间部分的 / pixlecount操作就应该可以变为 *Invpixlecount，其中Invpixlecount = 1.0f/pixlecount,变除法为乘法，而且这部分计算还可以很容易的用SSE实现。我测试过，改进后的实现和解析opencv中Box Filter的实现并提出进一步加速的方案（源码共享） 这篇文好章的速度不相上下，但这里有个优势就是可以并行的。另外，最重要的一点时，当要计算上述浮点版半径版本的boxfilter时，积分图是不需要再次重新计算的，而积分图的计算所占的耗时至少有一半左右。因此，这个场合使用积分图版本的盒子滤波会更有优势。
