上一篇文章談及了GIMP裡實現的小波分解，但是這僅僅是把影像分解為多層的資料，如果快速的獲取分解資料以及後續怎麼利用這些資料，則是本文的重點。
一、我們先來看看演算法速度的最佳化問題。
原始的GIMP實現需要將影像資料轉換為浮點數後，然後進行各級的模糊和圖層混合，這樣得到的結果是比較精確的，但是存在兩個方面的問題，一個是佔用了較多的記憶體，因為GIMP這個版本的小波分解各層是沒有改變資料的尺寸的，因此，如果使用浮點，佔用的記憶體要比位元組版本的大四倍，而且和層數有著密切的關係。第二個是浮點的處理還是稍微慢了點，雖然對現在的CPU來說，浮點數更易用SIMD指令集最佳化。但是如果有更好的資料型別的話，使用SIMD可以獲得更高的計算速度。
我們知道，位元組版本的模糊可以獲得很高的計算效率，這個場景的模糊是否可以使用呢，我個人分析認為，這個版本的演算法已經不適合用位元組版本的模糊了，原因主要是精度太低了，因為他每次的結果都和上一次的模糊相關，而我們透過GIMP的可是化介面可以看到，中間的每一個層很多畫素都是靠近127的，這就說明細節方面的資訊都是很小的數值。
個人覺得，對於這個演算法，我們可以把資料放大到unsigned short範圍，比如把原始的畫素值都擴大256倍，然後進行模糊，這樣模糊的精度就會大幅的提高，比如原始的9的畫素的加權累加值（歸一化後的權重）如果是100.3，那麼由於位元組版本取整的原因，最後的返回值為100，放大256倍後，則累加值變為25676.8，四捨五入取整後則為25677，而非100放大256倍的25600了，這樣多次模糊的精度就可以加以充分的保證（和浮點數比較還是有差異，但不影響結果）。
我們在稍微觀察下上一篇文章我們所提到的3*3的權重：
如果我們每個係數都放大16倍，我們會得到非常最佳化的一組權重係數：
可以看到，權重係數裡全是1、2、4，這種係數如果被乘數是整形的話，都可以直接用移位實現，而放大16倍最後的除法，也可以直接由右移實現，那這種計算過程就完全避免了乘法，效率可想而知可以得到極大的提高。
我們在考慮一種特別的最佳化，因為權重整形化後的累加值是16，那麼如果每個元素的最大值不超過4096，則累計值也就不會超過65536，這個時候如果是用普通C語言實現，其實沒有啥區別，但是我們知道SIMD指令確有所不同，他針對不同的資料型別有不同的乘法和加法指令，其所能計算的覆蓋範圍也不同，如果能用16位表達，則儘量用16位表達。因此，考慮影像資料的特殊性，如果我們只把他們的資料範圍擴大16倍，則不會超出4096的範圍，就可以滿足前面的這個假設。
放大16倍是否能滿足精度的需求呢，沒有具體理論的分析啊，個人覺得啊，至少在層數不大於4層時，差異不會很大，大於4層，也應該可以接受，留待實踐檢測吧。
對於上一篇文章所說到的取樣點位置的問題，我們必須考慮邊緣位置處的資訊，因為隨著取樣範圍的擴大，會有更多的取樣點超出影像的有效範圍，這個時候一個簡單的辦法就是實現準備好一副擴充套件過邊界的影像，這裡的擴充套件方法通道都選擇邊緣映象，而非重複邊緣。考慮到層數不會太多，擴充套件部分的邊緣計算量和佔用記憶體也不會太誇張。
同時，我們在考慮到演算法的特殊性，雖然取樣範圍很廣，但是真正用到的取樣點也只有9個，所以我們也可以使用類似於我在SSE影像演算法最佳化系列九：靈活運用SIMD指令16倍提升Sobel邊緣檢測的速度（4000*3000的24點陣圖像時間由480ms降低到30ms） 一文中使用到的技術，分配三行臨時的記憶體，每次計算前填充好三行的資料，不過注意一點，由於取樣的三行的記憶體並不是在行方向上連續的，因此，也不能像那個文章那樣，可以重複利用兩行的記憶體了，而必須每次都填充。
理論上說，這種填充技術要比前面的分配一整片的記憶體的記憶體填充工作量大3倍左右，速度應該要慢一些，優點是佔用的中間記憶體少一些，但是實測，還是這種的速度快一些，或許是因為三行記憶體的大小有效，訪問時的cache miss要好很多吧。
一個簡答的程式碼片段如下所示：
for (int Y = 0; Y < Height; Y++)
{
short *LinePD = Dest + Y * Width * Channel;
int PosY0 = ColOffset[Y] * Width * Channel, PosY1 = Y * Width * Channel, PosY2 = ColOffset[Y + Radius + Radius] * Width * Channel;
for (int X = 0, Index = 0; X < Radius; X++)
{
int PosX = RowOffset[X] * Channel;
memcpy(First + Index, Src + PosY0 + PosX, Channel * sizeof(unsigned short));
memcpy(Second + Index, Src + PosY1 + PosX, Channel * sizeof(unsigned short));
memcpy(Third + Index, Src + PosY2 + PosX, Channel * sizeof(unsigned short));
Index += Channel;
}
memcpy(First + Radius * Channel, Src + PosY0, Width * Channel * sizeof(unsigned short));
memcpy(Second + Radius * Channel, Src + PosY1, Width * Channel * sizeof(unsigned short));
memcpy(Third + Radius * Channel, Src + PosY2, Width * Channel * sizeof(unsigned short));
for (int X = 0, Index = (Width + Radius) * Channel; X < Radius; X++)
{
int PosX = RowOffset[X + Width + Radius] * Channel;
memcpy(First + Index, Src + PosY0 + PosX, Channel * sizeof(unsigned short));
memcpy(Second + Index, Src + PosY1 + PosX, Channel * sizeof(unsigned short));
memcpy(Third + Index, Src + PosY2 + PosX, Channel * sizeof(unsigned short));
Index += Channel;
}
int BlockSize = 8, Block = (Width * Channel) / BlockSize;
for (int X = 0; X < Block * BlockSize; X += BlockSize)
{
__m128i P0 = _mm_loadu_si128((__m128i *)(First + X));
__m128i P1 = _mm_loadu_si128((__m128i *)(First + X + Radius * Channel));
__m128i P2 = _mm_loadu_si128((__m128i *)(First + X + 2 * Radius * Channel));
__m128i P3 = _mm_loadu_si128((__m128i *)(Second + X));
__m128i P4 = _mm_loadu_si128((__m128i *)(Second + X + Radius * Channel));
__m128i P5 = _mm_loadu_si128((__m128i *)(Second + X + 2 * Radius * Channel));;
__m128i P6 = _mm_loadu_si128((__m128i *)(Third + X));
__m128i P7 = _mm_loadu_si128((__m128i *)(Third + X + Radius * Channel));
__m128i P8 = _mm_loadu_si128((__m128i *)(Third + X + 2 * Radius * Channel));
__m128i Sum0 = _mm_adds_epu16(_mm_adds_epu8(P0, P2), _mm_adds_epu16(P6, P8));
__m128i Sum1 = _mm_adds_epu16(_mm_adds_epu8(P1, P7), _mm_adds_epu16(P3, P5));
__m128i Sum2 = _mm_slli_epi16(P4, 2);
__m128i Sum = _mm_adds_epu16(_mm_adds_epu16(Sum0, Sum2), _mm_adds_epu16(Sum1, Sum1));
_mm_storeu_si128((__m128i *)(LinePD + X), _mm_srli_epi16(Sum, 4));
}
for (int X = Block * BlockSize; X < Width * Channel; X++)
{
LinePD[X] = (First[X] + First[X + (Radius + Radius) * Channel] + Third[X] + Third[X + (Radius + Radius) * Channel]
+ (First[X + Radius * Channel] + Second[X] + Second[X + (Radius + Radius) * Channel] + Third[X + (Radius) * Channel]) * 2
+ Second[X + Radius * Channel] * 4) / 16;
}
//for (int X = 0; X < Width; X++)
//{
// // 1 2 1
// // 2 4 2
// // 1 2 1
// for (int Z = 0; Z < Channel; Z++)
// {
// LinePD[Z] = (First[X * Channel + Z] + First[(X + Radius + Radius) * Channel + Z] + Third[X * Channel + Z] + Third[(X + Radius + Radius) * Channel + Z]
// + (First[(X + Radius) * Channel + Z] + Second[X * Channel + Z] + Second[(X + Radius + Radius) * Channel + Z] + Third[(X + Radius) * Channel + Z]) * 2
// + Second[(X + Radius) * Channel + Z] * 4) / 16;
// }
// LinePD += Channel;
//}
}
在小波分解的計算中，核心耗時的也就是這個模糊，其他的諸如圖層模式的混合，直接用SIMD指令也很是簡單。這裡不予以贅述。
二、小波資料的幾個簡單應用
（一）降噪
我們在搜尋GIMP的wavelet_decompse相關資訊時，搜到了一個叫krita的軟體，在其官網發現他有一個小波去噪的功能，而且這個軟體還是開源的，詳見網站https://krita.org/zh/。
稍微分析下其krita的程式碼，可以發現其分解的過程其實就是借鑑了GIMP的過程：
//main loop for(int level = 0; level < scales; ++level){ //copy original KisPaintDeviceSP blur = new KisPaintDevice(*original); //blur it KisWaveletKernel::applyWavelet(blur, rc, 1 << level, 1 << level, flags, 0); //do grain extract blur from original KisPainter painter(original); painter.setCompositeOpId(op); painter.bitBlt(0, 0, blur, 0, 0, rc.width(), rc.height()); painter.end(); //original is new scale and blur is new original results << original; original = blur; updater->setProgress((level * 100) / scales); }
精髓的1 << level也存在於這裡。
在其kis_wavelet_noise_reduction.cpp我們也可以看到這樣的程式碼：
for (float* it = begin; it < fin; it++) { if (*it > threshold) { *it -= threshold; } else if (*it < -threshold) { *it += threshold; } else { *it = 0.; } }
這裡有一個閾值的引數，即為外界客戶需要提供的引數。
後續我們在翻閱小波去噪的論文時，也多次發現類似的公式，其實這就是所謂的軟閾值處理：
那這裡的核心其實就是對小波分解後的每層資料，按其值大小進行一定的裁剪，注意這個裁剪最好不要處理Residual層。
我們嚴格的按照這個流程對GIMP的小波分解後的資料進行了降噪測試，其中幾個效果如下所示：
原圖 軟閾值，層數5，噪音閾值3 硬閾值，層數5，噪音閾值3
同樣引數下，軟閾值的去噪程度要比硬閾值大很多。
我們對照網上提供的matlab版本程式碼，似乎結果還是有蠻大的差異的。
close all; [A,map]=imread('C:\2.jpg'); x=rgb2gray(A); imshow(x); [c,s]=wavedec2(x,2,'sym4'); a1=wrcoef2('a',c,s,'sym4'); figure; imshow(uint8(a1)); a2=wrcoef2('a',c,s,'sym4',2); figure; imshow(uint8(a2));
我隨意拿了幾張人臉的圖去測試，結果意外發現，這個去噪和磨皮的效果有那麼幾分相似：
應該說，好好的把握處理好這幾個層的資料，應該還能有更多的結果出來，而且處理的自由度也比較高。
在處理速度上，預設5層資訊，3000*2000的灰度圖可以在90ms內處理完成，速度還是相當的快的。
(二)、銳化
上述去噪的過程中，我們將小於閾值的小波分量全部賦值為0了，理論上的目的是消除了噪音，將絕對值大於閾值的部分也進行了削弱，相當於減少了細節的資訊，那麼如果我們把這個過程稍微修改下，就可以產生很好的銳化效果。
我們這樣操作，設定兩個引數，一個Threshold，一個Amount引數，當小波係數小於Threshold，我們不做任何的處理，保留這部分，如果我們認為他是噪音，則表示噪音部分不做任何變動，如果大於Threshold，我們則放大小波係數，放大的程度取決於Amount引數， 這樣即增強了影像的細節部分，起到了銳化的作用，同時也不會過分的放大噪音，因此，就會比普通的銳化具有更強的識別性。
一種簡單的程式碼如下所示：
for (int X = 0; X < Width * Channel; X++) { if (LinePS[X] > Threshold) LinePS[X] += (LinePS[X] - Threshold) * Amount * 0.01; else if (LinePS[X] < -Threshold) LinePS[X] -= (Threshold - LinePS[X]) * Amount * 0.01; else { //LinePS[X] = 0; } LinePD[X] += LinePS[X]; }
由此產生的效果如下圖：
處理後的圖要銳利清晰很多。
（三）擴充套件
我們注意到這個小波的分解過程其實和我們的拉普拉斯金字塔的建立過程是非常類似的，只不過在拉普拉斯金字塔裡使用的5*5的一個加權模糊。而拉普拉斯金字塔裡儲存的也是類似於模糊之間數值的差異。因此上述的過程也可以直接適用於拉普拉斯金字塔處理。
使用拉普拉斯金字塔處理還有一個優勢就是速度可以進一步加快，畢竟其一層的尺寸是逐步變小的，處理量也就相對應的小了一些，比如處理3000*2000的灰度圖，5層金字塔耗時大概也就35ms。
三、小結
無論是小波分解，還是拉普拉斯分解，其更為重要的特點都是多尺度，那麼也可以將很多其他的單尺度的演算法放到這裡來，也會會有更多的意想不到的效果，特別是如果每一層的細節處理使用不同的自適應引數，可能會有更為廣闊的空間。
目前，我已經將小波去噪和小波銳化整合到我的SIMD最佳化的DEMO，詳見Enhance -> Denoise或者Enhance->Sharpen選單下。
有興趣的朋友可以從：https://files.cnblogs.com/files/Imageshop/SSE_Optimization_Demo.rar?t=1660121429 下載。
如果想時刻關注本人的最新文章，也可關注公眾號：
