拉普拉斯金字塔在多图HDR算法中的应用以及多曝光图像的融合算法简介。

这方面最经典的文章是2007年Tom Mertens等人发表的《Exposure Fusion》一文，用简单的篇幅和公式描述了一个非常优异的合成过程，虽然在2019年Charles Hessel发表了一篇《Extended Exposure Fusion》的文章中，提出了比Exposure Fusion更为优异的合成效果，但是代价是更高昂的计算成本，而Exposure Fusion也已经相当优秀了，本文主要简单记录下个人的Exposure Fusion优化过程。

　　Exposure Fusion的思路也非常之简单，输入是一系列图像对齐后的大小格式相同的图像，输出是一张合成的多细节图。那么在进行计算之前，他需要做以下的准备。 

　　1、对每副图像按照某些原则计算每融合的权重。文章里提出了3种权重。

　　   （1） 对比度：

               Contrast: we apply a Laplacian filter to the grayscale version of each image, and take the absolute value of the filter response 16. This yields a simple indicator C for contrast. It tends to assign a high weight to important elements such as edges and texture. 

　　对应的matlab代码非常简单：

% contrast measure
function C = contrast(I)
h = [0 1 0; 1 -4 1; 0 1 0]; % laplacian filter
N = size(I,4);
C = zeros(size(I,1),size(I,2),N);
for i = 1:N
    mono = rgb2gray(I(:,:,:,i));
    C(:,:,i) = abs(imfilter(mono,h,'replicate'));
end

　　我们可以认为就是个边缘检测。

　　（2）饱和度 

           Saturation: As a photograph undergoes a longer exposure, the resulting colors become desaturated and eventually clipped. Saturated colors are desirable and make the image look vivid. We include a saturation measure S, which is computed as the standard deviation within the R, G and B channel, at each pixel。

% saturation measure
function C = saturation(I)
N = size(I,4);
C = zeros(size(I,1),size(I,2),N);
for i = 1:N
    % saturation is computed as the standard deviation of the color channels
    R = I(:,:,1,i);
    G = I(:,:,2,i);
    B = I(:,:,3,i);
    mu = (R + G + B)/3;
    C(:,:,i) = sqrt(((R - mu).^2 + (G - mu).^2 + (B - mu).^2)/3);
end

　　也是非常简单的过程。

　　（3）曝光度

　　　Well-exposedness: Looking at just the raw intensities within a channel, reveals how well a pixel is exposed. We want to keep intensities that are not near zero (underexposed) or one (overexposed). We weight each intensity i based on how close it is to 0.5 using a Gauss curve: exp - (i-0.5)*(i-0.5)/(2*σ *σ), where σ equals 0.2 in our implementation. To account for multiple color channels, we apply the Gauss curve to each channel separately, and multiply the results, yielding the measure E。

% well-exposedness measure
function C = well_exposedness(I)
sig = .2;
N = size(I,4);
C = zeros(size(I,1),size(I,2),N);
for i = 1:N
    R = exp(-.5*(I(:,:,1,i) - .5).^2/sig.^2);
    G = exp(-.5*(I(:,:,2,i) - .5).^2/sig.^2);
    B = exp(-.5*(I(:,:,3,i) - .5).^2/sig.^2);
    C(:,:,i) = R.*G.*B;
end

　　每副图像得到三个指标（对比度C、饱和度S以及曝光度E），将他们相乘得到这幅图像的综合权重W。

     2、根据每副图像的权重，计算在序列中图像的每副图像的归一化权重，原文表述如下：

          To obtain a consistent result, we normalizethe values of the N weight maps such that they sum to one at each pixel (i, j): 

　　　注意，这是单个像素进行的处理，也就是说对于一个序列的图像，每个对应像素位置的权重先计算好后，再累加得到总权重，然后在每个对应像素的权重除以总权重得到归一化的值。

    3、理论上讲，得到了这些权重，就可以对N个图像进行直接融合，即使用下述公式：

 　　但是如果真的这样做，得到的结果惨不忍睹，即使我们对归一化后的权重进行高斯模糊、保边模糊等等也是解决不了问题的。所以这个文章的最核心的精华部分就在下面的过程中了。 

     4、我们并不是直接用这些权重进行合成，而是对原序列图像进行高斯拉普拉斯金字塔分解（利用其中的拉普拉斯数据），对各序列对应的归一化权重做高斯金字塔分解。 然后重构拉普拉斯金字塔，重构的规则为：

　　　　即对合成后各层拉普拉斯金字塔数据，用每个序列的权重的高斯金字塔数据乘以对应的拉普拉斯金字塔数据，然后累加。就如此简单。

　　关于这个算法的原理性说明，CSDN这个作者讲的也比较好： 【HDR】曝光融合（Exposure Fusion），有兴趣可以参考。

　　这个算法的效果通过测试确实还是可以的，原始的作者提供了MATLAB版本的代码，但是实际测试速度是非常慢的，这个也不怪这些论文的作者，他们的重点是实现算法本身，而不是工程化。

　　为了能提高算法的实用性，我们使用C++结合SIMD指令对该过程进行优化。优化的方式主要有以下几个方面：

       1、通过适当的改变数据类型来提高速度。

        M版本的代码是全程是使用的浮点类型的来计算各种数据的，为了效率，有些过程我们可以使用整形来代替。比如获取各个特征的权重的时，我们可以用byte来记录就可以了。在比如金字塔的分解和重构时，高斯金字塔我们也用byte类型记录，拉普拉斯金字塔考虑到有负数以及数据的范围，可以用signed short类型来记录。 

       2、通过改变数据类型后，使用对应的SIMD指令可以起到更大的步长，浮点数一般一次性只能处理4个数，如果是signed short则可以同时处理8个数据了。 

       3、对于金字塔分解和重构，要从原理上进行优化，详见：SSE图像算法优化系列二十六:和时间赛跑之优化高斯金字塔建立的计算过程。

       我们以计算饱和的过程为例，看看SSE的优化步骤： 

 　   简简单单的M代码，需要大概10倍以上行数的C代码和SIMD指令进行优化。 

　　实际上，我们在优化时，归一化的那个权重矩阵我也用byte类型来保存数据，原本以为这样每个图的权重有可能会比较小，存在精度损失，不过实际测试，和M代码的结果也没啥特别大的视觉区别。 

　　另外，还有内存方面的优化问题，如果建立所有图像的金字塔序列，然后再计算特征合成，这样会占用比较多的内存，特别是图像序列比较多时，实际上我们可以边分解边进行计算，这样带来的好处时速度有适当加速（应该还是cache miss较小的原因），当然如果每个序列都分配金字塔的内存，有一个好处就是，可以使用omp进行并行，适当的通过占用资源提高速度。

　　原始论文的作者在获取了各个特征的权重后，是使用的乘法来统计特征，我这里也使用了加法进行测试，测试结果是加法的结果比乘法的要稍微暗一点。区别也不是很大，但是要注意如果是使用乘法，有的时候有些特征的权重为0，为了不让其他特征的权重被淹没掉，建议加一后再乘。 

　　效果：  

　　其中第一个图是使用4副752*500的图合成，M代码大约耗时1600ms，我这里优化后的速度大概是25ms，有64倍的提速。 

       第二个图使用的是9副1024*683的图合成，M代码大约耗时5000ms，我这里优化后的速度大概是80ms，也有约60倍的提速。 

      从两幅图的合成结果来看，至少视觉上是非常不错的。不过在第一幅的中间靠上的部分的天空，以及第二副的右下角的台灯的中间部分还是有过曝的现象，这两个部分在原图中实际上是可以找到比较好的底图的，但是合成结果并没有吸收他们，这个可以适当的通过改变金字塔的数量来调节，但是没有啥理论依据去支持他，比如我们把他们的金字塔分别调整为6和7时效果如下两图所示：

　　这个现象在有些论文里也有提到，即金字塔层数太深，可能导致结果出现out-of-range的瑕疵，金字塔太少了，会有low-frequency halos出现。

       这个算法还有一个应用场景，在2019年的另外一篇论文里有提及： Simulated Exposure Fusion， 他的输入只有一张图片，然后通过某些手段生成一个序列的图像，然后在把这序列的图像予以使用EF算法合成，从而得到不错的结果，这个待我有空了也来研究一下，这也是个不错的主意。 

　　这个算法也应该非常适合使用GPU进行优化，那可能会获得更为客观的提速比，对于某些应用从而得到实时级别的效果。

　　在SSE图像算法优化系列二十九：基础的拉普拉斯金字塔融合用于改善图像增强中易出现的过增强问题（一）一文中使用的融合方法，实际上也是可以应用于多图的融合的，只不过这个时候低频的融合方式就不能是选择哪一个图了，此时可以使用平均值、最大值、最大值等策略，而高频融合一般只能使用AbsMax策略， 那个3*3 AbsMax方法也不可以使用了（因为一致性检测哪里无法实现）。

      不过这种融合对于多图来说一般效果不是很好，但是不排除也有其应用场景，他的一个优势就是速度能够做到更快，通常是本文放大的1/3耗时左右，而且内存的占用也很少。

　　另外，对于这个算法，还有一个现象就是其对图像的合成结果和用光度立体法对那些图的合成结果似乎很像，当然，仅仅是合成结果类似，广地立体能获得的其他信息这个算法是无法获取的。 

       提供一个测试DEMO供有兴趣的朋友玩玩：https://files.cnblogs.com/files/Imageshop/Exposure_Fusion.rar?t=1694501148&download=true

       或者有兴趣的朋友也可以在ipol网站中找到一个在线的DEMO以及代码，详见地址：https://ipolcore.ipol.im/demo/clientApp/demo.html?id=278