SSE图像算法优化系列二十九：基础的拉普拉斯金字塔融合用于改善图像加强中易出现的过加强问题（一）

　　拉普拉斯金字塔融合是多图融合相关算法里最简单和最容易实现的一种，咱们在看网络上大部分的文章都是在拿那个苹果和橙子融合在一块儿，变成一个果橙的效果做为例子说明。在这方面确实融合的比较好。可是本文咱们主要讲下这个在图像加强方面的运用。

       首先咱们仍是来说下这个融合的过程和算法优化。

       算法第一步：输入两个相同大小，位深的图像，经过拉普拉斯分解获得各自的拉普拉斯金字塔数据A和B。

　　算法第二步：选择下低频部分的融合规则，这里的低频部分，其实就是高斯金字塔最顶层那里的数据，这个数据至关因而原图像的一个高斯模糊的下采样版本，反应了基本的图像轮廓和信息。

       一般状况下，融合规则有三种：

       （1）选择A；

       （2）选择B；

       （3）选择A和B的平均值。

        算法第三步：选择高频部分的融合规则。高频表明了图像的边缘和细节，固然也多是噪音。高频部分的数据保存在各自拉普拉斯金字塔的除最顶层外的层中（最顶层和高斯金字塔的最顶层共享数据）。

       这里的融合规则就有多种，经常使用的好比以下几种：

       （1）选择A和B中绝对值最大的；

       （2）选择A和B领域中绝对值最大的。

　　第一种规则比较容易理解，绝对值大（拉普拉斯金字塔数据有多是负值得），表示这里的边缘强度越高，细节越丰富。

       第二种也好理解，一般咱们选择3*3领域。A的3*3领域的绝对值最大值若是大于B的3*3领域最大值，则选择A，不然选择B，这种作法的道理就是用领域去除必定的噪音影响。一般伴随着该种方法的还有一个叫作一致性检测的过程，即若是中心位置的融合系数选自原图像A变换的系数，而其周围领域内的融合系数大部分都选取自原图像B变换的系数 ，则把中心位置的融合系数修改成图像B变换后的系数，反之亦然。

   　那还有一种基于基于3X3窗口内类似性测度，获取拉普拉斯金字塔融合结果的方法，这个能够参考： https://wenku.baidu.com/view/c8ae11adf61fb7360b4c65c4.html ，这种融合规则因为考虑了与相邻像素间的相关性，下降了对边缘的敏感性，可以有效减小融合像素的错误选取，在必定程度上显著提升了融合算法的鲁棒性，从而提升了融合效果。 实测这种也还能够，可是代码比较麻烦，这里不描述。

　　算法第四步：高频和低频都已经处理好后，则重构图像获得结果值。

　   一个简单的描述过程以下：

int IM_LaplacePyramidFusion(unsigned char *SrcA, unsigned char *SrcB, unsigned char *Dest, int Width, int Height, int Stride, LowFrequencyFusionRule Low, HighFrequencyFusionRule High, int Level) { int Channel = Stride / Width; if ((SrcA == NULL) || (SrcB == NULL) || (Dest == NULL))                return IM_STATUS_NULLREFRENCE; if ((Width <= 0) || (Height <= 0))                                    return IM_STATUS_INVALIDPARAMETER; if ((Channel != 1) && (Channel != 3) && (Channel != 4))                return IM_STATUS_INVALIDPARAMETER; int Status = IM_STATUS_OK; Level = IM_ClampI(Level, 1, IM_GetMaxPyramidLevel(Width, Height, 1));        // 通过测试若是直接处理到最小为1个像素的金子塔，效果并很差
 Pryamid *GaussPyramid = (Pryamid *)calloc(Level, sizeof(Pryamid));            // 必须用calloc，否则在后面的释放函数中可能存在野指针释放问题，高斯金字塔能够用同一个内存
    Pryamid *LaplacePyramidA = (Pryamid *)calloc(Level, sizeof(Pryamid));        // 图A的拉普拉斯金字塔
    Pryamid *LaplacePyramidB = (Pryamid *)calloc(Level, sizeof(Pryamid));        // 图B的拉普拉斯金字塔
    Pryamid *GaussPyramidD, *LaplacePyramidD; if ((GaussPyramid == NULL) || (LaplacePyramidA == NULL) || (LaplacePyramidB == NULL)) { Status = IM_STATUS_OUTOFMEMORY; goto FreeMemory; } Status = IM_AllocatePyramidMemory(GaussPyramid, Width, Height, Stride, Level, true, sizeof(unsigned char));            // 分配内存，高斯金字塔的塔底就是原数据
    if (Status != IM_STATUS_OK)    goto FreeMemory; Status = IM_AllocatePyramidMemory(LaplacePyramidA, Width, Height, Stride, Level, false, sizeof(unsigned char)); if (Status != IM_STATUS_OK)    goto FreeMemory; Status = IM_AllocatePyramidMemory(LaplacePyramidB, Width, Height, Stride, Level, false, sizeof(unsigned char)); if (Status != IM_STATUS_OK)    goto FreeMemory; GaussPyramid[0].Data = SrcA; for (int Y = 1; Y < Level; Y++)        // 各级高斯金字塔
 { Status = IM_DownSample8U((unsigned char *)GaussPyramid[Y - 1].Data, (unsigned char *)GaussPyramid[Y].Data, GaussPyramid[Y - 1].Width, GaussPyramid[Y - 1].Height, GaussPyramid[Y - 1].Stride, GaussPyramid[Y].Width, GaussPyramid[Y].Height, GaussPyramid[Y].Stride, Channel); if (Status != IM_STATUS_OK)  goto FreeMemory; } // 拉普拉斯金子塔的最顶层和高斯金字塔的是同样的额
    memcpy(LaplacePyramidA[Level - 1].Data, GaussPyramid[Level - 1].Data, LaplacePyramidA[Level - 1].Height * LaplacePyramidA[Level - 1].Stride); for (int Y = Level - 2; Y >= 0; Y--) { Status = IM_UpSampleSub8U((unsigned char *)GaussPyramid[Y + 1].Data, (unsigned char *)GaussPyramid[Y].Data, (unsigned char *)LaplacePyramidA[Y].Data, GaussPyramid[Y + 1].Width, GaussPyramid[Y + 1].Height, GaussPyramid[Y + 1].Stride, GaussPyramid[Y].Width, GaussPyramid[Y].Height, GaussPyramid[Y].Stride, Channel); if (Status != IM_STATUS_OK)  goto FreeMemory; } GaussPyramid[0].Data = SrcB; for (int Y = 1; Y < Level; Y++)        // 各级高斯金字塔
 { Status = IM_DownSample8U((unsigned char *)GaussPyramid[Y - 1].Data, (unsigned char *)GaussPyramid[Y].Data, GaussPyramid[Y - 1].Width, GaussPyramid[Y - 1].Height, GaussPyramid[Y - 1].Stride, GaussPyramid[Y].Width, GaussPyramid[Y].Height, GaussPyramid[Y].Stride, Channel); if (Status != IM_STATUS_OK)  goto FreeMemory; } // 拉普拉斯金子塔的最顶层和高斯金字塔的是同样的额
    memcpy(LaplacePyramidB[Level - 1].Data, GaussPyramid[Level - 1].Data, LaplacePyramidB[Level - 1].Height * LaplacePyramidB[Level - 1].Stride); for (int Y = Level - 2; Y >= 0; Y--) { Status = IM_UpSampleSub8U((unsigned char *)GaussPyramid[Y + 1].Data, (unsigned char *)GaussPyramid[Y].Data, (unsigned char *)LaplacePyramidB[Y].Data, GaussPyramid[Y + 1].Width, GaussPyramid[Y + 1].Height, GaussPyramid[Y + 1].Stride, GaussPyramid[Y].Width, GaussPyramid[Y].Height, GaussPyramid[Y].Stride, Channel); if (Status != IM_STATUS_OK)  goto FreeMemory; } LaplacePyramidD = LaplacePyramidA; // 低频部分的融合
    PyramidLowFreqFusion((unsigned char *)LaplacePyramidA[Level - 1].Data, (unsigned char *)LaplacePyramidB[Level - 1].Data, (unsigned char *)LaplacePyramidD[Level - 1].Data, LaplacePyramidA[Level - 1].Width, LaplacePyramidA[Level - 1].Height, LaplacePyramidA[Level - 1].Stride, Low); for (int Y = 0; Y < Level - 1; Y++)                    // 高频部分的融合
 { if (High == SinglePixelAbsMax) PyramidHighFreqFusion_AbsMax((unsigned char *)LaplacePyramidA[Y].Data, (unsigned char *)LaplacePyramidB[Y].Data, (unsigned char *)LaplacePyramidD[Y].Data, LaplacePyramidA[Y].Width, LaplacePyramidA[Y].Height, LaplacePyramidA[Y].Stride); else if (High == LocalAbsMaxWithConsistencyCheck) PyramidHighFreqFusion_3X3MaxAbsValue((unsigned char *)LaplacePyramidA[Y].Data, (unsigned char *)LaplacePyramidB[Y].Data, (unsigned char *)LaplacePyramidD[Y].Data, LaplacePyramidA[Y].Width, LaplacePyramidA[Y].Height, LaplacePyramidA[Y].Stride); // else //PyramidHighFreqFusion_3X3Similarity(LaplacePyramidA[Y], LaplacePyramidB[Y], LaplacePyramidD[Y], PryamidW[Y], PryamidH[Y], Channel);*/
 } GaussPyramidD = GaussPyramid; memcpy(GaussPyramidD[Level - 1].Data, LaplacePyramidD[Level - 1].Data, LaplacePyramidD[Level - 1].Height * LaplacePyramidD[Level - 1].Stride); for (int Y = Level - 2; Y >= 0; Y--)            // 重构拉普拉斯金子塔
 { if (Y != 0) IM_UpSampleAdd8U((unsigned char *)GaussPyramidD[Y + 1].Data, (unsigned char *)LaplacePyramidD[Y].Data, (unsigned char *)GaussPyramidD[Y].Data, GaussPyramidD[Y + 1].Width, GaussPyramidD[Y + 1].Height, GaussPyramidD[Y + 1].Stride, GaussPyramidD[Y].Width, GaussPyramidD[Y].Height, GaussPyramidD[Y].Stride, Channel); else IM_UpSampleAdd8U((unsigned char *)GaussPyramidD[Y + 1].Data, (unsigned char *)LaplacePyramidD[Y].Data, Dest, GaussPyramidD[Y + 1].Width, GaussPyramidD[Y + 1].Height, GaussPyramidD[Y + 1].Stride, GaussPyramidD[Y].Width, GaussPyramidD[Y].Height, GaussPyramidD[Y].Stride, Channel); } FreeMemory: IM_FreeGaussPyramid(GaussPyramid, Level, true); IM_FreeLaplacePyramid(LaplacePyramidA, Level); IM_FreeLaplacePyramid(LaplacePyramidB, Level); if (GaussPyramid != NULL) free(GaussPyramid); if (LaplacePyramidA != NULL) free(LaplacePyramidA); if (LaplacePyramidB != NULL) free(LaplacePyramidB); return Status; }

　　咱们上面的全部的高斯或者拉普拉斯金字塔数据都是用unsigned char类型来描述的， 为何能够这样作呢，作个简单的分析。第一，高斯金字塔用byte是毫无疑问的，第二，前面说了，严格的拉普拉斯金字塔是有负数的，可是咱们考虑到一个这个负数大于-127的可能性是很是小的，这种状况可能会在二值图像中出现，而二值图的处理算法中能用到金字塔吗，我彷佛没据说过，因此咱们能够把拉普拉斯金字塔的数据加上127，让总体在0和255之间，这样有不少算法都直接调用了。

       在咱们的高频或者低频的选取过程当中，由于都不存在新的数据出来，也就是没有啥几何乘积计算，所以，用byte保存也不存在啥大问题。

　　PyramidLowFreqFusion低频部分的融合代码很是简单，以下所示:

/// 低频部分的融合，通常有三种方式。一、取图像A的系数； 二、取图像B的系数；三、取图像A和个B系数的平均值。 /// 其实这里应该用高斯金字塔的最顶层数据，只是因为拉普拉斯和高斯金字塔共享这一层数据，因此也能够直接这样写
int PyramidLowFreqFusion(unsigned char *LaplacePyramidA, unsigned char *LaplacePyramidB, unsigned char *LaplacePyramidD, int Width, int Height, int Stride, LowFrequencyFusionRule Low) { int Channel = Stride / Width; if ((LaplacePyramidA == NULL) || (LaplacePyramidB == NULL) || (LaplacePyramidD == NULL))    return IM_STATUS_NULLREFRENCE; if ((Width <= 0) || (Height <= 0))                                                            return IM_STATUS_INVALIDPARAMETER; if ((Channel != 1) && (Channel != 3) && (Channel != 4))                                        return IM_STATUS_INVALIDPARAMETER; if (Low == SrcA) memcpy(LaplacePyramidD, LaplacePyramidA, Height * Stride); else if (Low == SrcB) memcpy(LaplacePyramidD, LaplacePyramidB, Height * Stride); else { // 也能够考虑某一副图占比高一点 //int BlockSize = 16, Block = (Height * Stride) / BlockSize; //for (int Y = 0; Y < Block * BlockSize; Y += BlockSize) //{ // __m128i SrcA = _mm_loadu_si128((__m128i *)(LaplacePyramidA + Y)); // __m128i SrcB = _mm_loadu_si128((__m128i *)(LaplacePyramidB + Y)); // __m128i Dst1 = _mm_srli_epi16(_mm_add_epi16(_mm_mullo_epi16(_mm_cvtepu8_epi16(SrcA), _mm_set1_epi16(3)), _mm_cvtepu8_epi16(SrcA)), 2); // __m128i Dst2 = _mm_srli_epi16(_mm_add_epi16(_mm_mullo_epi16(_mm_cvtepu8_epi16(_mm_srli_si128(SrcA, 8)), _mm_set1_epi16(3)), _mm_cvtepu8_epi16(_mm_srli_si128(SrcA, 8))), 2); // _mm_storeu_si128((__m128i *)(LaplacePyramidD + Y), _mm_packus_epi16(Dst1, Dst2)); //} //for (int Y = Block * BlockSize; Y < Height * Stride; Y++) //{ // LaplacePyramidD[Y] = (LaplacePyramidA[Y] * 3 + LaplacePyramidB[Y]) >> 2; // 最高层(低频）系数取平均 //} // **************************** 真正意义上的平均值 *************************************
        int BlockSize = 16, Block = (Height * Stride) / BlockSize; for (int Y = 0; Y < Block * BlockSize; Y += BlockSize) { _mm_storeu_si128((__m128i *)(LaplacePyramidD + Y), _mm_avg_epu8(_mm_loadu_si128((__m128i *)(LaplacePyramidA + Y)), _mm_loadu_si128((__m128i *)(LaplacePyramidB + Y)))); } for (int Y = BlockSize * BlockSize; Y < Height * Stride; Y++) { LaplacePyramidD[Y] = (LaplacePyramidA[Y] + LaplacePyramidB[Y]) >> 1;    // 最高层(低频）系数取平均
 } } return IM_STATUS_OK; }

　　取平均直接用_mm_avg_epu8就能够了。

　　高频部分若是选择绝对值最大的方案代码也是很简单的：

/// 基于系数绝对值取大的融合策略进行拉普拉斯金字塔图像融合。
int PyramidHighFreqFusion_AbsMax(unsigned char *LaplacePyramidA, unsigned char *LaplacePyramidB, unsigned char *LaplacePyramidD, int Width, int Height, int Stride) { int Channel = Stride / Width; if ((LaplacePyramidA == NULL) || (LaplacePyramidB == NULL) || (LaplacePyramidD == NULL))    return IM_STATUS_NULLREFRENCE; if ((Width <= 0) || (Height <= 0))                                                            return IM_STATUS_INVALIDPARAMETER; if ((Channel != 1) && (Channel != 3) && (Channel != 4))                                        return IM_STATUS_INVALIDPARAMETER; int BlockSize = 16, Block = (Height * Stride) / BlockSize; __m128i C127 = _mm_set1_epi8(127); for (int Y = 0; Y < Block * BlockSize; Y += BlockSize) { __m128i SrcA = _mm_loadu_si128((__m128i *)(LaplacePyramidA + Y)); __m128i SrcB = _mm_loadu_si128((__m128i *)(LaplacePyramidB + Y)); __m128i Flag = _mm_cmpgt_epu8(_mm_absdiff_epu8(SrcA, C127), _mm_absdiff_epu8(SrcB, C127)); _mm_storeu_si128((__m128i *)(LaplacePyramidD + Y), _mm_blendv_epi8(SrcB, SrcA, Flag)); } for (int Y = Block * BlockSize; Y < Height * Stride; Y++) { if (IM_Abs(LaplacePyramidA[Y] - 127) > IM_Abs(LaplacePyramidB[Y] - 127)) LaplacePyramidD[Y] = LaplacePyramidA[Y]; else LaplacePyramidD[Y] = LaplacePyramidB[Y]; } return IM_STATUS_OK; }

　　其中的_mm_absdiff_epu8函数以下所示：

// 返回8位字节数数据的差的绝对值
inline __m128i _mm_absdiff_epu8(__m128i a, __m128i b) { return _mm_or_si128(_mm_subs_epu8(a, b), _mm_subs_epu8(b, a)); }

　　这里要减去127的主要缘由是前面所说的再计算拉普拉斯金字塔时咱们增长了127，而这里计算时咱们须要真正的拉普拉斯金子塔数据。

       用_mm_blendv_epi8能够方便的解决后续的抉择问题，有点至关于C语言的里的三目运算符。

　　关于PyramidHighFreqFusion_3X3MaxAbsValue这个函数就要复杂不少了，首先这种3*3的领域计算，我仍是推荐我在sobel边缘算子优化一文中提到的那种优化结构，能够支持In-Place操做，同时还能够完美处理边缘。算法的流程是标准化的。就是求出各自3*3领域的绝对值的最大值，而后进行比较，为了后续的一致性检测，比较的结果须要写入个临时内存，在实现时，咱们作了以下处理：

_mm_storeu_si128((__m128i *)(LinePF + X), _mm_cmpgt_epu8(MaxA, MaxB));

　　其中的MaxA和MaxB为领域的最大值，这里也就是说A>B，对应的Flag位置设置为255，不然设置为0（也是用的unsigned char内存保存的）。

　　为何这样作，有两个好处，第一，咱们在后续的一致性检测里，能够充分利用这个数据的特殊性。在一致性检测里，咱们要判断周边的是否是大部分都和中心点来自同一个数据源，一种处理方式就是把周边的8个点的数据都相加，若是中心点为0，周边的和大于255*4，则代表周边和中心不太一致，须要把中心的点改成255，若是中心点为255，而周边的点的和小于255*4，则中心点要改成0，用代码表示以下：

int Sum = First[X + 0] + First[X + 1] + First[X + 2] + Second[X + 0] + Second[X + 2] + Third[X + 0] + Third[X + 1] + Third[X + 2];
if (LinePD[X] == 0 && Sum > 255 * 4)            //    若是当前点为0，而且周边8个点中至少有5个点为1，则把当前点的值修改成1。
    LinePD[X] = 255;
else if (LinePD[X] == 255 && Sum < 255 * 4)     //    若是当前点为1，而且周边8个点中至少有5个点为0，则把当前点的值修改成0。
    LinePD[X] = 0;

    　　第二，在SSE优化时，这个特殊性是能够帮上大忙的，主要体如今两个方面，一时如上的Sum计算过程，咱们若是直接按照255相加，则8个数会超出8位所表达的范围，这样就要转换到16位的空间进行计算了，可是若是咱们把epu8当作epi8，这个时候255就编程了-1，此时的加法我使用_mm_add_epi8，则结果能在epi8的范围内，这样一次性就能够处理16个像素了。二是后续咱们须要根据这个Flag对组中的输出结果作明示，这个时候咱们就能够直接使用这个Flag作mask供_mm_blendv_epi8调用。

　　咱们在俩看看上面的判断部分如何用SSE处理，由于SSE不善于作分支，因此咱们须要想办法，这样作，咱们看看下面的代码是否是和上面的一个意思：

//    if ((LinePD[X] == 0 && Sum > 255 * 4) || ((LinePD[X] == 255 && Sum < 255 * 4)))
//        LinePD[X] = 255 - LinePD[X];

 　　可是这里是有不一样的，咱们能够很方便的对上述代码SSE处理：

__m128i FlagA = _mm_and_si128(_mm_cmpeq_epi8(Current, _mm_setzero_si128()), _mm_cmplt_epi8(Sum, _mm_set1_epi8(-4)));
__m128i FlagB = _mm_and_si128(_mm_cmpeq_epi8(Current, _mm_set1_epi8(255)), _mm_cmpgt_epi8(Sum, _mm_set1_epi8(-4)));
__m128i FlagAB = _mm_or_si128(FlagA, FlagB);
_mm_storeu_si128((__m128i *)(LinePD + X), _mm_blendv_epi8(Current, _mm_subs_epu8(_mm_set1_epi8(255), Current), FlagAB));    //    局部取反

　　这样效率能够大大的提升。

　　优化方面基本讲完了，固然，因为我没有共享代码，大部分实际上是写给本身看的，由于我怕时间长了本身都不知道为何要这样写。

       下面咱们测试下算法效果和性能：

    　　　　　　　　　　　　　　　　　　SrcA　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　SrcB

  

                                             低频=SrcA，高频=3*3领域, Level = 10　　　　　　　　　　　　　　　　　　　　　　　　 低频=SrcB，高频=3*3领域, Level = 10

  

 　　        　 低频=(SrcA + SrcB) / 2，高频=3*3领域, Level = 10　　　　　　　　　　　　　　　　　　　　 低频=(SrcA + SrcB) / 2，高频=3*3领域, Level =5

　　上述原图B是某个加强算法处理后的结果，很明显，改算法对图像右下角的暗部的加强效果很好，可是同时图像上部的天空区域已经彻底过曝了，天空的云消失不见了，这在不少加强算法中都会出现相似的状况，而在原图中天空的细节自己就已经比较好了，所以，咱们尝试用不一样的选项对这两幅图作拉普拉斯融合，若是高频选择SrcA，则总体融合后的图像暗部加强的不明显，选择SrcB，则天空恢复的不够好，选择(SrcA + SrcB) / 2则能对天空和暗部都有较好的恢复。

        另外，金字塔的层数对结果也有必定的影响，在最后两张图中，咱们能够看到Level=5时的效果要比为10时的稍微好一点，咱们通常也不建议高斯金字塔的最顶层取得过小，一般，咱们取5层金字塔应该能得到较为满意的效果。

　　效率方面，通常1920*1080的彩色图像，这种混合大概须要20-30ms左右（取决于选择的参数），一半的时间用于了金字塔的构建。

       其余说明：

       一、在PyramidLowFreqFusion函数中，咱们注释掉了一部分，这部分注释的代码的本意是低频的算法咱们不必定必定要是取平均值，也能够根据实际的状况更增强调某一对象，好比假如SrcB是处理后的部分，咱们能够把他的权重设置为75%，而SrcA的权重设置为25%。

       二、对于彩色图像，若是三个通道独立写，则对每一个像素，有可能每一个通道的高频或低频部分会选自不一样的来源，这样有可能致使结果出现异常的彩色，一种解决方案是采用高频或低频部分的灰度信号做为判断源。

      三、金子塔融合的基本原理仍是保留更多的细节，所以，若是对一个原始图像进行了相似锐化方面的处理后，这个图在和原图进行融合，那基本上不会有什么变化的，柔和的结果必然是靠近锐化后的结果图。这个你们能够本身作实验。

      四、这个融合的过程基本不须要外接的参数接入，咱们能够考虑把他做为某些算法的最后一个默认步骤。

      五、对于任意两幅大小相同的图，这个算法融合的结果也是蛮有意思的，以下：

  

    　固然，这种融合仍是有必定的限制的，下一节，咱们将讨论基于蒙版的金字塔融合，那里能够更加智能的获取更好的融合效果。

        提供一个DEMO供测试效果：极度优化版本工程：https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar，见MultiImage->LaplacePyramidFusion菜单。