SSE图像算法优化系列二十八：深度优化局部拉普拉斯金字塔滤波器。

基于局部拉普拉斯金字塔的Edge-aware滤波器是在2011年由Adobe 公司的研究员Sylvain Paris（大神级人物，写了很多文章）提出的，我在4年前曾经参考有关代码实现过这个算法，但是速度也是非常慢的，所以当时也没有继续做深入的研究，前段时间做另外一个算法时仔细的研究了下高斯和拉普拉斯金子塔的优化，因此又抽时间仔细的分析了算法的论文和代码，由于论文的理论部分还有一些我没有想清楚，因此在这里我只对研读过程中涉及的代码方面的优化做个解读。

　　经过我最终的优化，处理1920 * 1024的彩色图，在保证效果不会有明显的瑕疵的取样值的情况下，大概能获得60ms的速度。

　　先分享下参考资料：

　　（1）论文 Local Laplacian Filters: Edge-aware Image Processing with a Laplacian Pyramid。　

　　（2）论文 Fast Local Laplacian Filters: Theory and Applications

　　（3）函数：matlab 2017的locallapfilt函数。

　　（4）插件：Paint.net的Laplacian pyramid filter effect （可反编译为C#代码)

　我们先看下源作者给出的一个最原始的matlab的过程：

 1 % naive O(N ^ 2) version for reference
 2 
 3 G = gaussian_pyramid(I);                                        %   由输入图像I建立高斯金字塔序列
 4 L = laplacian_pyramid(zeros(size(I)));                             %   建立一个和图像I大小相匹配的空的拉普拉斯金字塔
 5 for lev0 = 1:length(L) - 1                                        %    对每一层金字塔
 6     for y0 = 1 : size(G{ lev0 }, 1)                                %    遍历每一层金字塔的高度
 7         for x0 = 1 : size(G{ lev0 }, 2)                            %    遍历每一层金字塔的宽度
 8             g0 = G{ lev0 }(y0, x0, :);                            %   每一层高斯金字塔的每个数据
 9             Iremap = r(I, g0);                                    %   根据每个数据值，计算输入图像I对应的一个映射新图像
10             Lremap = laplacian_pyramid(Iremap, lev0 + 1);        %   计算这个新映射图像的拉普拉斯金子塔
11             L{ lev0 }(y0, x0, :) = Lremap{ lev0 }(y0, x0, :);   %   填空前面建立的空的拉普拉斯金字塔对应位置的值
12         end
13     end
14 end
15 L{ end } = G{ end };
16 R = reconstruct_laplacian_pyramid(L);                            %   根据拉普拉斯金字塔构建出结果图像

　　原文核心的关于这段代码的解释如下（为了不浪费空间，下面的图片是从原文截图，然后用PS再把短行拼接成长行的）：

　　简单的说，就是需要遍历所有高斯金字塔图像中的所有像素，根据每个像素的像素值，都由原图和某个映射函数重新计算出一个和原图一样大小的图像，然后计算该图像的拉普拉斯金字塔，如上述代码的第10行所示，注意此时的拉普拉斯金字塔只需要构建到对应的像素所在的高斯金字塔那一层就可以了，然后呢，取该像素位置所对应临时金字塔的值作为结果图像在此位置的拉普拉斯金字塔值。当所有层的像素都计算完成后，结果图的拉普拉斯金子塔就构建完成了，这个时候结果图像就可以由拉普拉斯金字塔重构出。

　　由上述分析可见，直接实现这个过程将是非常耗时的，每一层金字塔的每一个系数都要靠构造一次拉普拉斯金字塔，如果图像宽和高都为N，则理论上说，所有的金字塔加起来是有N*N+N*N个像素的，这个时候就需要计算N*N大小图像的拉普拉斯金字塔（N*N+N*N)次，会让算法根本无法使用。在原始论文中，作者提到为了计算某个位置的拉普拉斯金字塔值，并不需要计算整体的值，而只需要取某个局部区域来计算，可以将计算的复杂度降低到O(N * Log(N))，确实如此，但是这样的过程依旧很慢很慢。

　　在没有看Fast Local Laplacian Filters: Theory and Applications论文之前，我想到的关于此方法的优化手段非常有效，因为对于常规8位图像来说，其像素的可能值只有0到255之间的256个值，因此，在上述N*N+N*N次构建拉普拉斯金字塔的过程中，最多其实只会有256种不同结果，这也就意味着我们只需要构建256次拉普拉斯金字塔就可以得到所有的结果。大概的matlab代码如下所示：

G = gaussian_pyramid(I);                                        
L = laplacian_pyramid(zeros(size(I)));                             
for M = 0 : 255
    Iremap = r(I, M);                                    
    Lremap = laplacian_pyramid(Iremap, length(L) + 1);
    for lev0 = 1:length(L) - 1                                        
        for y0 = 1 : size(G{ lev0 }, 1)                                
            for x0 = 1 : size(G{ lev0 }, 2)                            
                if (G{ lev0 }(y0, x0, :) = M)                    
                    L{ lev0 }(y0, x0, :) = Lremap{ lev0 }(y0, x0, :) 
                end 
            end 
        end
    end
end
L{ end } = G{ end };
R = reconstruct_laplacian_pyramid(L);

　　结合后面将要介绍的处理方法，利用C++和SSE处理一幅1920 * 1024（2M Pixel)的灰度图，单线程大概的耗时在1200ms左右，而源作者的论文使用8核并行计算处理1M像素的图都用了4000ms多，提速了很多倍，处理效果则是一摸一样。

　　而论文Fast Local Laplacian Filters: Theory and Applications则做了更多的做工，他首先分析局部拉普拉斯算法和双边滤波的关系，然后分析了这个算法慢的主要原因，最后提出了他自己的解决方案，正如上面我们的所分析的，我们只需要做256次完整的拉普拉斯分解就可以了，而根据采样定理，其实不一定要做这么多次，只要多于某个采样数值时，系统一样可以稳定的输出，而这个数值通常都要远远的小于256次，当然，我们是不太方便直接从图像中找到这个最小的取样数据的，但是经过摸索，一般来说大于12次以上效果都是不错的，这篇论文的作者提出的相关参考代码如下：

 1 % INPUT
 2 % I : input greyscale image
 3 % r : a function handle to the remaping function
 4 % N : number discretisation values of the intensity
 5 % OUTPUT
 6 % F : filtered image
 7 function [F]=llf(I,sigma,fact,N)
 8     [height width]=size(I);
 9     n_levels=ceil(log(min(height,width))-log(2))+2;
10     discretisation=linspace(0,1,N);
11     discretisation_step=discretisation(2);
12     input_gaussian_pyr=gaussian_pyramid(I,n_levels);
13     output_laplace_pyr=laplacian_pyramid(I,n_levels);
14     output_laplace_pyr{n_levels}=input_gaussian_pyr{n_levels};
15     for ref=discretisation
16         I_remap=fact*(I-ref).*exp(-(I-ref).*(I-ref)./(2*sigma*sigma));
17         temp_laplace_pyr=laplacian_pyramid(I_remap,n_levels);
18         for level=1:n_levels-1
19             output_laplace_pyr{level}=output_laplace_pyr{level}+...
20                 (abs(input_gaussian_pyr{level}-ref)<discretisation_step).*...
21                 temp_laplace_pyr{level}.*...
22                 (1-abs(input_gaussian_pyr{level}-ref)/discretisation_step);            
23         end
24     end
25     F=reconstruct_laplacian_pyramid(output_laplace_pyr);

　　其中第19到22行即为插值的过程，我们测试了一下对应的Matlab代码，处理1M像素的图，大概耗时在3800ms左右，我们知道matlab的代码确实只适合教学，因此，优化的余地是有很大的。

　在优化前，我们还是定性的说下上面过程中涉及到的reampping Function，在原始的论文中，作者提到了这个函数起到了细节和边缘调整的作用，对于高斯金字塔中的任一像素值g₀，我们设定一个参数б_r , 当原图I中的像素i的值在g₀附近时，我们认为这些点属于g₀附近的细节，而远离g₀的部分则属于边缘，对细节和边缘我们采用两个不同的处理函数r_d和r_e，一般要求r_d和r_e必须是单调递增函数，而且满足，也及时要求r_d和r_e在分界处是连续的。

　　一个典型的处理如下所示：

else

　　其中在做细节处理时，论文建议的f_d和f_e如下所示：

其对应的曲线如下所示：

简单的分析下图片的直观认识吧，我们看看detail smoothing的曲线，在输入为g₀时输出为g_0，在小于g₀的б_r 范围内，输出是大于输入的，而在大于g₀的б_r范围内，输出是小于于输入的。也就是说在g₀附近的像素是朝向g₀进一步靠近的，从而使得这一块的细节都趋于一致，而在远离g₀的位置，像素未受到影响，这样在整体的表现上如果g₀在原始图像中属于某个平坦区域，则其周边的像素也慢慢往g₀靠近，如果他属于边缘，则周边的像素基本未改变，这个时候通过金字塔则可以把这种改变通过领域的关系一步一步的代入到拉普拉斯系数中，这也是所谓的局部带来的好处。在detail enhancement的增强中，则是一个相反的过程。

但在Fast Local Laplacian Filters: Theory and Applications一文配套的代码中，我们看到作者并没有使用上述曲线来处理，而是使用了这样的一个函数：

其中sigma一般取值0.15左右，f为用户调节参数，但f小于0时，起到平滑作用，大于0时，起到细节增强作用。

　　注意上述公式中的i和g₀都必须是归一化后的值。

　　下左图示出了此曲线（高斯曲线）和原始论文曲线（指数曲线）的差异，其中g₀取值为0.5。可见此曲线在整个定义域是非常光滑的，在远离g₀处并没有像原始论文那样对像素毫无影响，但影响也是非常小了。

　　上图的右侧部分为f取不同值时的曲线分布情况，我们注意到当f=-2或者f=4时的曲线出现了异常情况，他已经不符合原始论文提出的曲线是单调递增函数的要求，此时图像也会出现一些异常情况，后续会给出这个异常的结果图。

　　那么实际上我们还有很多其他的曲线可以使用，比如有关代码里列出如下的曲线函数：

%This is just a toy example!
function y=remapping_function(x)
   % y=(x-0.1).*(x>0.1)+(x+0.1).*(x<-0.1); %smoothing
    y=3.*x.*(abs(x)<0.1)+(x+0.2).*(x>0.1)+(x-0.2).*(x<-0.1); %enhancement
end

　　实际测试也是没有问题的，也能达到类似的效果。

　好了，下面我们来集中力量来实现上述新代码的C++优化。首先，为了较为准确的实现这个过程，我们先把图像数据转换为浮点数。但是我们可能不做归一化的处理，即浮点的范围还是控制在0和255之间。那么金字塔建立的优化过程再次不在赘述，可以参考我的相关博客。核心的优化就在第16到22行代码，我们先来处理第16行代码：

　　　　I_remap=fact*(I-ref).*exp(-(I-ref).*(I-ref)./(2*sigma*sigma))；

　　我们将他翻译为普通的C++代码：

float Inv = 1.0f / 255.0f / 255.0f / (2 * sigma * sigma);;
for (int I = 0; I < GaussPyramid[0].Height * GaussPyramid[0].Width; I++)
{
    float I = GaussPyramid[0][I], Diff = I - ref;
    GaussPyramidT[I] = IM_ClampF(I + f * Diff * exp(-Diff * Diff * Inv), 0, 255);
}

  GaussPyramid[0]其实就是原始输入图像，GaussPyramidT是临时的高斯金字塔数据，IM_ClampF是个裁剪限幅函数，主要是为了让数据不产生溢出，实际测试好像不限幅也没什么大的问题。
  这段代码是非常耗时的，第一，他是对原始图大小进行处理的，第二，exp的计算是个非常缓慢的过程，因此，我们尝试了SSE处理。

int BlockSize = 4, Block = (GaussPyramid[0].Height * GaussPyramid[0].Width) / BlockSize;
for (int I = 0; I < Block * BlockSize; I += BlockSize)
{
    __m128 V = _mm_loadu_ps(GaussPyramid[0] + I);
    __m128 Diff = _mm_sub_ps(V, _mm_set1_ps(ref));
    __m128 Exp = _mm_fexp_ps(_mm_neg_ps(_mm_mul_ps(_mm_mul_ps(Diff, Diff), _mm_set1_ps(Inv))));
    __m128 Dst = _mm_add_ps(_mm_mul_ps(_mm_mul_ps(Diff, Exp), _mm_set1_ps(f)), V);
    Dst = _mm_min_ps(_mm_max_ps(Dst, _mm_setzero_ps()), _mm_set1_ps(255.0f));
    _mm_storeu_ps(GaussPyramidT[0] + I, Dst);
}
for (int I = Block * BlockSize; I < GaussPyramid[0].Height * GaussPyramid[0].Width; I++)
{
    float I = GaussPyramid[0][I], Diff = I - ref;
    GaussPyramidT[I] = IM_ClampF(I + f * Diff * exp(-Diff * Diff * Inv), 0, 255);
}

　　基本上是直接翻译，其中有些函数为自定义的，如下所示：

//    求负数
inline __m128 _mm_neg_ps(__m128 a)
{
    return _mm_sub_ps(_mm_setzero_ps(), a);
}
//    http://martin.ankerl.com/2007/10/04/optimized-pow-approximation-for-java-and-c-c/
//    快速的指数运算，精度一般
inline __m128 _mm_fexp_ps(__m128 x)
{
    __m128i T = _mm_cvtps_epi32(_mm_add_ps(_mm_mul_ps(x, _mm_set1_ps(1512775)), _mm_set1_ps(1072632447)));
    __m128i TL = _mm_unpacklo_epi32(_mm_setzero_si128(), T);
    __m128i TH = _mm_unpackhi_epi32(_mm_setzero_si128(), T);
    return _mm_movelh_ps(_mm_cvtpd_ps(_mm_castsi128_pd(TL)), _mm_cvtpd_ps(_mm_castsi128_pd(TH)));
}

　　我们测试发现，虽然使用的是精度一般版本的exp函数，但是对最终的结果影响并不大，但是速度的提升则非常明显，而所谓的Clamp处理则可以直接使用min和max函数实现。

　　接着就是第18到22行代码，我们直接翻译为普通C语言如下所示：

for (level = 1; level < n_levels - 1; level++)
{
    for (int I = 0; I < (PryamidW[level] * PryamidH[level]; I++)
    {
        if (abs(GaussPyramid[level][I] - ref) < discretisation_step)            //    插值计算，表示在有效的范围内
        {
            LaplacePyramid[level][I] += (LaplacePyramidT[level][I] * (1 - abs(GaussPyramid[level][I] - ref) / discretisation_step));
        }
    }
}

　　重点是处理内部的循环。

　　因为内部循环里有条件判断和跳转，一般来说是不利于SSE处理的，如果要用SSE处理，其实施要点是找到某个参数，是的跳转和不跳转通过该参数能用同一个公式计算，我们观察上式，但不符合那个判断跳转时，我们可以通过把LaplacePyramidT[level][I]这个变量用0值代替，这样可保证结果不变。另外一个可以考虑的地方就是，如果存在较多的多数据同时不满足条件的情况，可以使用_mm_movemask_ps的函数来做处理，如果他返回值为0，我们可以不继续后续的处理，否则，就统一处理，如下所示：

int BlockSize = 4, Block = (PryamidW[level] * PryamidH[level]) / BlockSize;
float *D = LaplacePyramid[level];
float *S = GaussPyramid[level];
float *T = LaplacePyramidT[level];
for (int I = 0; I < Block * BlockSize; I += BlockSize)
{
    __m128 Abs = _mm_abs_ps(_mm_sub_ps(_mm_loadu_ps(S + I), _mm_set1_ps(ref)));
    __m128 Cmp = _mm_cmplt_ps(Abs, _mm_set1_ps(Step));
    if (_mm_movemask_ps(Cmp) != 0)
    {
        __m128 LT = _mm_blendv_ps(_mm_setzero_ps(), _mm_loadu_ps(T + I), Cmp);
        _mm_storeu_ps(D + I, _mm_sub_ps(_mm_add_ps(_mm_loadu_ps(D + I), LT), _mm_mul_ps(LT, _mm_mul_ps(Abs, _mm_set1_ps(1.0f / discretisation_step)))));
    }
}
for (int I = Block * BlockSize; I < PryamidW[level] * PryamidH[level]; I++)
{
    // do something
}

　　其中用到的一些自定义函数如下:

//    浮点数据的绝对值计算
inline __m128 _mm_abs_ps(__m128 v)
{
    static const int i = 0x7fffffff;
    float mask = *(float*)&i;
    return _mm_and_ps(v, _mm_set_ps1(mask));
}

　　通过以上优化方式，我们最终的测试结果表明针对2M大小的图像，平均处理时间稳定在145ms左右（取样数设置为12），相对原先的情况已经有了非常大的提高。

　　在耗时组成方面，我们测试数据如下，临时的高斯-拉普拉斯金字塔的构建85ms，映射函数耗时30ms，填充拉普拉斯金字塔数据30ms。可见大部分时间还是用在金字塔的处理上。

　　我们已经知道，整数版本的（8位或者32位的）金字塔的建立要比浮点版本快很多，特备是8位的金字塔数据，如果我们使用整数版本的效果会如何呢，我们进行了使用8位金子塔的版本进行了测试，当然，为了精度，拉普拉斯金字塔的数据部分还是需要使用singed short类型来保存，测试的效果表明，整数版本的精度也是足够的。如下图所示：

　　　　　　　　　　　　　　原图　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　浮点版本　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　8位金字塔版本

　　仔细对比，会有一点点的差异，但是基本上靠眼睛是区分不出来的。

　　但是，使用8位金字塔的第一个优势是建立金字塔的速度非常快，第二，函数的映射部分，我们可以使用查找表的方式快速实现（取整数结果），第三，填充拉普拉斯金字塔这一块使用相关整数处理的SSE指令，也可以一次性处理8个像素了，因此都大为加速，综合这三个优势，我们最终的优化结果做到了2M像素60ms的处理速度。

　　在Fast Local Laplacian Filters: Theory and Applications一文中，作者也给出了他优化的处理速度，如下所示：

　　可见，同比，我的速度比他的实现快了很多很多倍，当然还是没有赶上GPU的速度，但是也相差不大，比如我的速度折算到4M像素，需要120ms，他的GPU版本比如的快了2.5倍左右，但是我用的是单线程的，如果考虑多线程，还是有一定的提速空间（虽然笨算法不易多线程了）。

　　再者，我们还是来讨论下映射函数的问题，前面讲了，快速版本的代码使用的映射函数并没有使用原始论文的版本，所以我们尝试把这个替换一下，得到的结论就是，原始版本的映射函数不适合插值使用，效果如下所示：

　　　　　　　　　　　　原图　　　　　　　　　　　　　　　　　　　　　　　　　采样数为12(б_r=0.16, alpha = 0.3)时的效果　　　　　　　　　　　　　　　采样为32(б_r=0.16, alpha = 0.3)时的效果

　　　　　　直接实现的效果（256色阶）　　　　　　　　　　　　　　　　　　　　　　　　　带防止噪音扩大的效果

　　可见这个映射函数在采样率较小时的效果是非常差的，具有明显的色块效果，而只有不进行任何下采样时效果才较为满意。

　　而注意到上面的第四个图，可以看到在原来平坦的区域里出现很多不起眼的噪点，在原始论文里作者也注意到了这个问题，他提出了一个解决方案，即限制最小的差异的放大程度，当我们需要增强细节时（alpha < 1），使用一个混合公式来修正函数f_d的值。

　　式中的系数T值由abs(i-g₀)/бr 决定，当该值小于0.01时，为0，当大于0.02时，为1，而介于两者之间是使用一个平滑函数修正，这样做的结果就是使得在和g₀特别接近时，相关的像素不会得到修正，而噪音的引起也就是因为这些特别相近的像素经过原来的处理后会变得差异很大引起的。所以这样做就可以有效地避免该问题，一个常用的平滑函数如下所示：

inline float IM_SmoothStep(const float t)
{
    return t * t * (3 - 2 * t);　　　　　　　　//　　powf(t, 2.0f) * powf(t - 2, 2.0f);
}

　此时的修正曲线如下图所示：

　　至于为什么原始论文的曲线不适合采样，我分析可能还是因为他的曲线是由两个函数组合而成的，而中间的曲线部分在采样时不能够获得足够的取样点，而造成数据丢失，而改进后的论文中的曲线是一条光滑连续的曲线，其曲率变化也非常自然，很少的取样点就可以获得较为理想的效果。同时，我们也主要到高斯曲线的映射结果似乎不存在噪音过增强的现象，而且也不需要采用类似指数映射那种处理方式来减少该问题，如果采用，反而可能会对结果图像带来光晕。

　　接下来我们分析另外一个问题，现在我们推荐使用高斯曲线来进行数据的映射，当函数中f取值小于0时，是处于一个去燥或者说平滑图像的作用，同时还能有效地保留边缘，当f大于0时，起到了细节增强或者说锐化的作用，因此，f小于0时该滤波器的作用相当于一个保边滤波器，比如他也可以有效地用在磨皮中（f = -1左右）：

　　　但是当f过分的小或者过分的大时，比如f=-2或f=4时，如上述曲线所示，此时的高斯曲线已经不是单调递增函数了，此时的图像会出现什么效果，我们做了测试：

　　　　　　　　原图　　　　　　　　　　　　　　　　　　　　　　　　　　　　f = -2　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　f = 4

　　我们看看中间这幅图，可以看到天空中的云原先的白色已经变成了灰色了，但是整个图像还是处于平滑的状态，明显处于一个结果错误的状态，而后面一个图的增强程度似乎很过，但是相对于中间的部分而言要稍微合理一点。

　　引起该结果的一个核心原因正如上述所言，映射函数出了问题，在此时的映射函数不同的两个输入，可以有相同的输出。所以在使用高斯曲线时一定要注意这个问题。

　　其他细节：

　　在论文里还提到其他的一些细节，比如说为了提高速度，我们可以只对底层的若干层拉普拉斯金字塔做上述操作，而更高层的金字塔数据不做处理，也可以只对高层的拉普拉斯金字塔进行重构，而底层的不做处理，特别是底层的不错处理，会极大的提高速度，毕竟第二层的数据只有第一层数据的1/4了，那效果如何了，我们下面给出了一些结果：

　　　　　　　　　　　　　　　　原图　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 f =2, 0层金字塔不处理　　　　　　　　　　　　　　　　　　　　　f = 4, 0层金字塔不处理

　　可见这种修改虽能加快速度，但效果不好，虽有一定的细节增强效果，但0层的影响太强烈。在看另外一个效果：

　　　　　　　　　　　　　　　　原图　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 0层不处理（f = -1）　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　0层处理（f = -1)

　　中间的结果是0层不处理的结果，我们惊喜的发现这个结果和简单探讨可牛影像软件中具有肤质保留功能的磨皮算法及其实现细节一文中的有几分的相似，我们分析认为当0层不处理时，原图的纹理就保存原始0层的拉普拉斯金字塔中，而0层以后的数据都进行了平滑的处理，这样数据在最后的叠合后就呈现了纹理保留（肤质保留）的功能，这也是所谓的巧合吧。

我们之前也描述几篇保边滤波器的文章，他们通常只具有保边效果，而这里的拉普拉斯滤波器确内在的可以实现保边和细节增强的作用，而且改变使用不同的映射函数还可以实现tone mapping、style transefer等较为高级的功能，不失为一个开创性的算法，我个人觉得他从传统的一些Edge-aware算法中能脱颖而出，利用最常见和简单的金字塔算法实现通用的效果，真的有点类似当然何凯明的暗通道算法一样。不过好像从作者发表该论文后，还少有作者进行进一步的算法拓展和应用，这也是比较郁闷之处啊。

　　最后，我们还拿这个算法和其他算法的效果做个对比，来看看这个算法的强大。