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算法中能脫穎而出，利用最多見和簡單的金字塔算法實現通用的效果，真的有點相似固然何凱明的暗通道算法同樣。不過好像從做者發表該論文後，還少有做者進行進一步的算法拓展和應用，這也是比較鬱悶之處啊。

　　最後，咱們還拿這個算法和其餘算法的效果作個對比，來看看這個算法的強大。

          

　　　　　　　　　　　　　　　　原圖　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　USM銳化（半徑=20，Amount =200, Threshold = 0）    

      

　　　　　　　　基於導向濾波的銳化（數量 = 200）　　　　　　　　　　　　　　　　　　　　　　　　　　局部金字塔濾波（採樣數=12，f = 2)    

　　能夠看到，局部金字塔在邊緣保留這一塊作的特別的好。

　　對於醫學圖像，該算法也能很好的起到加強做用：

      

　　原始的渾濁圖像通過處理後變得很是清晰。

　　醫學圖像一般有不少都是16位的，該算法對16位依然有效，只是處理過程要稍做修改便可。

　   Demo下載地址：https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar