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菜單。