matlab中表示拉普拉斯分布_深度优化局部拉普拉斯金字塔滤波器。
微信公眾號：OpenCV學堂關注獲取更多計算機視覺與深度學習知識
覺得文章有用，請戳底部【好看】支持
算法概述
基于局部拉普拉斯金字塔的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的過程：
?%?naive?O(N?^?2)?version?for?referenceG?=?gaussian_pyramid(I);????????????????????????????????????????%???由輸入圖像I建立高斯金字塔序列
L?=?laplacian_pyramid(zeros(size(I)));?????????????????????????????%???建立一個和圖像I大小相匹配的空的拉普拉斯金字塔
for?lev0?=?1:length(L)?-?1????????????????????????????????????????%????對每一層金字塔
????for?y0?=?1?:?size(G{?lev0?},?1)????????????????????????????????%????遍歷每一層金字塔的高度
????????for?x0?=?1?:?size(G{?lev0?},?2)????????????????????????????%????遍歷每一層金字塔的寬度
????????????g0?=?G{?lev0?}(y0,?x0,?:);????????????????????????????%???每一層高斯金字塔的每個數據
????????????Iremap?=?r(I,?g0);????????????????????????????????????%???根據每個數據值，計算輸入圖像I對應的一個映射新圖像
????????????Lremap?=?laplacian_pyramid(Iremap,?lev0?+?1);????????%???計算這個新映射圖像的拉普拉斯金子塔
????????????L{?lev0?}(y0,?x0,?:)?=?Lremap{?lev0?}(y0,?x0,?:);???%???填空前面建立的空的拉普拉斯金字塔對應位置的值
????????end
????end
end
L{?end?}?=?G{?end?};
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次以上效果都是不錯的，這篇論文的作者提出的相關參考代碼如下：
%?INPUT%?I?:?input?greyscale?image
%?r?:?a?function?handle?to?the?remaping?function
%?N?:?number?discretisation?values?of?the?intensity
%?OUTPUT
%?F?:?filtered?imagefunction?[F]=llf(I,sigma,fact,N)
????[height?width]=size(I);
????n_levels=ceil(log(min(height,width))-log(2))+2;
????discretisation=linspace(0,1,N);
????discretisation_step=discretisation(2);
????input_gaussian_pyr=gaussian_pyramid(I,n_levels);
????output_laplace_pyr=laplacian_pyramid(I,n_levels);
????output_laplace_pyr{n_levels}=input_gaussian_pyr{n_levels};
????for?ref=discretisation
????????I_remap=fact*(I-ref).*exp(-(I-ref).*(I-ref)./(2*sigma*sigma));
????????temp_laplace_pyr=laplacian_pyramid(I_remap,n_levels);
????????for?level=1:n_levels-1
????????????output_laplace_pyr{level}=output_laplace_pyr{level}+...
????????????????(abs(input_gaussian_pyr{level}-ref)????????????????temp_laplace_pyr{level}.*...
????????????????(1-abs(input_gaussian_pyr{level}-ref)/discretisation_step);????????????
????????end
????end
????F=reconstruct_laplacian_pyramid(output_laplace_pyr);
其中第19到22行即為插值的過程，我們測試了一下對應的Matlab代碼，處理1M像素的圖，大概耗時在3800ms左右，我們知道matlab的代碼確實只適合教學，因此，優化的余地是有很大的。
在優化前，我們還是定性的說下上面過程中涉及到的reampping Function，在原始的論文中，作者提到了這個函數起到了細節和邊緣調整的作用，對于高斯金字塔中的任一像素值g0，我們設定一個參數бr , 當原圖I中的像素i的值在g0附近時，我們認為這些點屬于g0附近的細節，而遠離g0的部分則屬于邊緣，對細節和邊緣我們采用兩個不同的處理函數rd和re，一般要求rd和re必須是單調遞增函數，而且滿足，也及時要求rd和re在分界處是連續的。一個典型的處理如下所示
簡單的分析下圖片的直觀認識吧，我們看看detail smoothing的曲線，在輸入為g0時輸出為g0，在小于g0的бr 范圍內，輸出是大于輸入的，而在大于g0的бr 范圍內，輸出是小于于輸入的。也就是說在g0附近的像素是朝向g0進一步靠近的，從而使得這一塊的細節都趨于一致，而在遠離g0的位置，像素未受到影響，這樣在整體的表現上如果g0在原始圖像中屬于某個平坦區域，則其周邊的像素也慢慢往g0靠近，如果他屬于邊緣，則周邊的像素基本未改變，這個時候通過金字塔則可以把這種改變通過領域的關系一步一步的代入到拉普拉斯系數中，這也是所謂的局部帶來的好處。在detail enhancement的增強中，則是一個相反的過程。
但在Fast Local Laplacian Filters: Theory and Applications一文配套的代碼中，我們看到作者并沒有使用上述曲線來處理，而是使用了這樣的一個函數：
其中sigma一般取值0.15左右，f為用戶調節參數，但f小于0時，起到平滑作用，大于0時，起到細節增強作用。注意上述公式中的i和g0都必須是歸一化后的值。
下左圖示出了此曲線(高斯曲線)和原始論文曲線(指數曲線)的差異，其中g0取值為0.5。可見此曲線在整個定義域是非常光滑的，在遠離g0處并沒有像原始論文那樣對像素毫無影響，但影響也是非常小了。
上圖的右側部分為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;?%smoothing
????y=3.*x.*(abs(x)<0.1)+(x+0.2).*(x>0.1)+(x-0.2).*(x0.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?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處理。
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?1;?level++){
????for?(int?I?=?0;?I?????{
????????if?(abs(GaussPyramid[level][I]?-?ref)?????????{
????????????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?{
????__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?{
????//?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位金字塔的第一個優勢是建立金字塔的速度非常快，第二，函數的映射部分，我們可以使用查找表的方式快速實現(取整數結果)，第三，填充拉普拉斯金字塔這一塊使用相關整數處理的SSE指令，也可以一次性處理8個像素了，因此都大為加速，綜合這三個優勢，我們最終的優化結果做到了2M像素60ms的處理速度。
在Fast Local Laplacian Filters: Theory and Applications一文中，作者也給出了他優化的處理速度，如下所示：
可見，同比，我的速度比他的實現快了很多很多倍，當然還是沒有趕上GPU的速度，但是也相差不大，比如我的速度折算到4M像素，需要120ms，他的GPU版本比如的快了2.5倍左右，但是我用的是單線程的，如果考慮多線程，還是有一定的提速空間(雖然笨算法不易多線程了)。
再者，我們還是來討論下映射函數的問題，前面講了，快速版本的代碼使用的映射函數并沒有使用原始論文的版本，所以我們嘗試把這個替換一下，得到的結論就是，原始版本的映射函數不適合插值使用，效果如下所示：
可見這個映射函數在采樣率較小時的效果是非常差的，具有明顯的色塊效果，而只有不進行任何下采樣時效果才較為滿意。
而注意到上面的第四個圖，可以看到在原來平坦的區域里出現很多不起眼的噪點，在原始論文里作者也注意到了這個問題，他提出了一個解決方案，即限制最小的差異的放大程度，當我們需要增強細節時(alpha < 1)，使用一個混合公式來修正函數fd的值。
式中的系數T值由abs(i-g0)/бr 決定，當該值小于0.01時，為0，當大于0.02時，為1，而介于兩者之間是使用一個平滑函數修正，這樣做的結果就是使得在和g0特別接近時，相關的像素不會得到修正，而噪音的引起也就是因為這些特別相近的像素經過原來的處理后會變得差異很大引起的。所以這樣做就可以有效地避免該問題，一個常用的平滑函數如下所示：
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時，如上述曲線所示，此時的高斯曲線已經不是單調遞增函數了，此時的圖像會出現什么效果，我們做了測試：
我們看看中間這幅圖，可以看到天空中的云原先的白色已經變成了灰色了，但是整個圖像還是處于平滑的狀態，明顯處于一個結果錯誤的狀態，而后面一個圖的增強程度似乎很過，但是相對于中間的部分而言要稍微合理一點。
引起該結果的一個核心原因正如上述所言，映射函數出了問題，在此時的映射函數不同的兩個輸入，可以有相同的輸出。所以在使用高斯曲線時一定要注意這個問題。
其他細節：
在論文里還提到其他的一些細節，比如說為了提高速度，我們可以只對底層的若干層拉普拉斯金字塔做上述操作，而更高層的金字塔數據不做處理，也可以只對高層的拉普拉斯金字塔進行重構，而底層的不做處理，特別是底層的不錯處理，會極大的提高速度，畢竟第二層的數據只有第一層數據的1/4了，那效果如何了，我們下面給出了一些結果：
可見這種修改雖能加快速度，但效果不好，雖有一定的細節增強效果，但0層的影響太強烈。在看另外一個效果：
中間的結果是0層不處理的結果，我們驚喜的發現這個結果和簡單探討可牛影像軟件中具有膚質保留功能的磨皮算法及其實現細節 一文中的有幾分的相似，我們分析認為當0層不處理時，原圖的紋理就保存原始0層的拉普拉斯金字塔中，而0層以后的數據都進行了平滑的處理，這樣數據在最后的疊合后就呈現了紋理保留(膚質保留)的功能，這也是所謂的巧合吧。
我們之前也描述幾篇保邊濾波器的文章，他們通常只具有保邊效果，而這里的拉普拉斯濾波器確內在的可以實現保邊和細節增強的作用，而且改變使用不同的映射函數還可以實現tone mapping、style transefer等較為高級的功能，不失為一個開創性的算法，我個人覺得他從傳統的一些Edge-aware算法中能脫穎而出，利用最常見和簡單的金字塔算法實現通用的效果，真的有點類似當然何凱明的暗通道算法一樣。不過好像從作者發表該論文后，還少有作者進行進一步的算法拓展和應用，這也是比較郁悶之處啊。
效果比對
最后，我們還拿這個算法和其他算法的效果做個對比，來看看這個算法的強大。
可以看到，局部金字塔在邊緣保留這一塊做的特別的好。
對于醫學圖像，該算法也能很好的起到增強作用：
原始的渾濁圖像經過處理后變得非常清晰。
醫學圖像通常有很多都是16位的，該算法對16位依然有效，只是處理過程要稍作修改即可。
Demo下載地址：https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar
作者網名：laviewpbt
是圖像處理，算法實現與加速優化方面的大神！其開發的imageshop軟件大小只有1MB，卻實現了非常豐富與復雜的各種圖像處理功能，他的郵箱地址為：Email: laviewpbt@sina.com
總結
以上是生活随笔為你收集整理的matlab中表示拉普拉斯分布_深度优化局部拉普拉斯金字塔滤波器。的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： 地平线机器人_地平线机器人CEO余凯：基
- 下一篇： iservice封装有哪些方法_5w大功