基于局部拉普拉斯金字塔的Edge-aware濾波器是在2011年由Adobe 公司的研究員Sylvain Paris（大神級人物，寫了很多文章）提出的，我在4年前曾經參考有關代碼實作過這個算法，但是速度也是非常慢的，是以當時也沒有繼續做深入的研究，前段時間做另外一個算法時仔細的研究了下高斯和拉普拉斯金子塔的優化，是以又抽時間仔細的分析了算法的論文和代碼，由于論文的理論部分還有一些我沒有想清楚，是以在這裡我隻對研讀過程中涉及的代碼方面的優化做個解讀。
基于局部拉普拉斯金字塔的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的過程：
原文核心的關于這段代碼的解釋如下（為了不浪費空間，下面的圖檔是從原文截圖，然後用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代碼如下所示：
結合後面将要介紹的處理方法，利用C++和SSE處理一幅1920 * 1024（2M Pixel)的灰階圖，單線程大概的耗時在1200ms左右，而源作者的論文使用8核并行計算處理1M像素的圖都用了4000ms多，提速了很多倍，處理效果則是一摸一樣。
而論文Fast Local Laplacian Filters: Theory and Applications則做了更多的做工，他首先分析局部拉普拉斯算法和雙邊濾波的關系，然後分析了這個算法慢的主要原因，最後提出了他自己的解決方案，正如上面我們的所分析的，我們隻需要做256次完整的拉普拉斯分解就可以了，而根據采樣定理，其實不一定要做這麼多次，隻要多于某個采樣數值時，系統一樣可以穩定的輸出，而這個數值通常都要遠遠的小于256次，當然，我們是不太友善直接從圖像中找到這個最小的取樣資料的，但是經過摸索，一般來說大于12次以上效果都是不錯的，這篇論文的作者提出的相關參考代碼如下：
其中第19到22行即為插值的過程，我們測試了一下對應的Matlab代碼，處理1M像素的圖，大概耗時在3800ms左右，我們知道matlab的代碼确實隻适合教學，是以，優化的餘地是有很大的。
在優化前，我們還是定性的說下上面過程中涉及到的reampping Function，在原始的論文中，作者提到了這個函數起到了細節和邊緣調整的作用，對于高斯金字塔中的任一像素值g0，我們設定一個參數бr , 當原圖I中的像素i的值在g0附近時，我們認為這些點屬于g0附近的細節，而遠離g0的部分則屬于邊緣，對細節和邊緣我們采用兩個不同的處理函數rd和re，一般要求rd和re必須是單調遞增函數，而且滿足
，也及時要求rd和re在分界處是連續的。
一個典型的處理如下所示：
else
其中在做細節處理時，論文建議的fd和fe如下所示：
其對應的曲線如下所示：
簡單的分析下圖檔的直覺認識吧，我們看看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時的曲線出現了異常情況，他已經不符合原始論文提出的曲線是單調遞增函數的要求，此時圖像也會出現一些異常情況，後續會給出這個異常的結果圖。
那麼實際上我們還有很多其他的曲線可以使用，比如有關代碼裡列出如下的曲線函數：
實際測試也是沒有問題的，也能達到類似的效果。
好了，下面我們來集中力量來實作上述新代碼的C++優化。首先，為了較為準确的實作這個過程，我們先把圖像資料轉換為浮點數。但是我們可能不做歸一化的處理，即浮點的範圍還是控制在0和255之間。那麼金字塔建立的優化過程再次不在贅述，可以參考我的相關部落格。核心的優化就在第16到22行代碼，我們先來處理第16行代碼：
I_remap=fact*(I-ref).*exp(-(I-ref).*(I-ref)./(2*sigma*sigma))；
我們将他翻譯為普通的C++代碼：
基本上是直接翻譯，其中有些函數為自定義的，如下所示：
我們測試發現，雖然使用的是精度一般版本的exp函數，但是對最終的結果影響并不大，但是速度的提升則非常明顯，而所謂的Clamp處理則可以直接使用min和max函數實作。
接着就是第18到22行代碼，我們直接翻譯為普通C語言如下所示：
重點是處理内部的循環。
因為内部循環裡有條件判斷和跳轉，一般來說是不利于SSE處理的，如果要用SSE處理，其實施要點是找到某個參數，是的跳轉和不跳轉通過該參數能用同一個公式計算，我們觀察上式，但不符合那個判斷跳轉時，我們可以通過把LaplacePyramidT[level][I]這個變量用0值代替，這樣可保證結果不變。另外一個可以考慮的地方就是，如果存在較多的多資料同時不滿足條件的情況，可以使用_mm_movemask_ps的函數來做處理，如果他傳回值為0，我們可以不繼續後續的處理，否則，就統一處理，如下所示：
其中用到的一些自定義函數如下:
通過以上優化方式，我們最終的測試結果表明針對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），使用一個混合公式來修正函數fd的值。
式中的系數T值由abs(i-g0)/бr 決定，當該值小于0.01時，為0，當大于0.02時，為1，而介于兩者之間是使用一個平滑函數修正，這樣做的結果就是使得在和g0特别接近時，相關的像素不會得到修正，而噪音的引起也就是因為這些特别相近的像素經過原來的處理後會變得差異很大引起的。是以這樣做就可以有效地避免該問題，一個常用的平滑函數如下所示：
此時的修正曲線如下圖所示：
至于為什麼原始論文的曲線不适合采樣，我分析可能還是因為他的曲線是由兩個函數組合而成的，而中間的曲線部分在采樣時不能夠獲得足夠的取樣點，而造成資料丢失，而改進後的論文中的曲線是一條光滑連續的曲線，其曲率變化也非常自然，很少的取樣點就可以獲得較為理想的效果。同時，我們也主要到高斯曲線的映射結果似乎不存在噪音過增強的現象，而且也不需要采用類似指數映射那種處理方式來減少該問題，如果采用，反而可能會對結果圖像帶來光暈。
接下來我們分析另外一個問題，現在我們推薦使用高斯曲線來進行資料的映射，當函數中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