【短道速滑十】非區域性均值濾波的指令集最佳化和加速（針對5*5的搜尋特例，可達到單核1080P灰度圖 28ms/幀的速度）。

Imageshop發表於2023-10-08

非區域性均值濾波(Non Local Means）作為三大最常提起來的去燥和濾波演算法之一（雙邊濾波、非區域性均值、BM3D），也是有著很多的論文作為研究和比較的物件，但是也是有著致命的缺點，速度慢，嚴重的影響了演算法的應用範圍。目前在已有的文獻中尚未看到在不對演算法的本質原理上做更改的情況下，能取得實時的效果，本文呢，也不求得到這個目的，只是對現有的開放的資源上來取得更進一步的提升。

　　標準的NL-Means演算法中，一般有三個引數，搜尋半徑SearchRadius，塊半徑PatchRadius，以及一個決定平滑程度的高斯函式引數Delta。在百度上能夠搜尋到的大部分文章所描述的提速演算法都是使用積分圖來提升NL-Means的速度，這也是目前來說唯一比較靠譜的最佳化技術，透過積分圖，可以做到演算法和塊半徑PatchRadius的大小基本無關，和Delta也無關，和SearchRadius成平方關係。

　　因此，在我們很多的嚴重的噪音影像中，SearchRadius至少需要取到7以上（涉及15*15= 225個領域範圍）才有明顯的效果，因此，這就相當於要計算225次全圖的某種計算，即使是每次這種計算只需要1ms（通常，任何影像處理演算法無法超越同等記憶體大小的memcpy的大小的），也需要225ms的，因此，確實比較感尷。

　　透過多執行緒方式可以適當對這個過程進行加速，畢竟每個畫素點的處理相對來說還是獨立的，但是，這個加速也收到物理核心的限制，就是8核的機器，滿利用，也無法達到8倍的加速效果。

話說回來，這麼大的計算量，用GPU也都是很吃力的。

　　目前使用積分圖來記性NL-Means演算法的比較好的文章也是大家常看的還是這一篇：

　　　　　非區域性均值濾波（NL-means）演算法的積分圖加速原理與C++實現

　　其實積分圖一直有個問題，可能很多搞影像的人都沒有注意到，或者說這個問題可能對某些演算法的影響還不是很大，不足有對大家注意到或者關注到，即積分圖存在著一下2個方面的問題和缺點：

　　1、當影像較大時，積分圖無法使用int型別來儲存，我們必須選擇能夠容納更大資料範圍的資料型別來儲存。

　　如果我們儲存的是一副位元組影像的積分圖，考慮極端情況，每個影像的值都是255，則積分影像最多隻可儲存uint.Maxvalue / 255 = 16843009個畫素，大約4100*4100大小的影像，一般來說，這個規模對於實際的應用來說是足夠了。

　　但是我們在實際中用到的很多情況，不是直接的影像積分圖，還常用的有平方積分圖，即求影像平方值後的積分圖，這個時候每一個元素的最大值達到了65025，極限情況下uint型別只可儲存uint.Maxvalue / 255 =66051個元素的總和，只大概是指256*256個影像的大小了，已經遠遠的無法滿足實際的應用需求。

　　因此，為了實現這個要求，我們必須選擇能夠容納更大資料範圍的資料型別，這裡有三個選擇：long long(int64) / float / double

　　第一個long long型別，即64位的整形資料，這個資料的表達範圍已經完全夠我們在影像中使用積分圖了，而且儲存的資料是非常準確的，特別是對於影像方面的資料來說（都是整形的），但是有個致命的問題：速度相當相當的慢，特別是同樣的計算和int型別比較的話，那真的不是一個檔次上的。我一直不太理解，現在大部分都是64位系統，為什麼對64位的資料的支援還是這麼的弱。而且我們看大部分指令集最佳化的函式對64位整形的支援都比較少。因此，非常不建議使用這個型別。

　　第二個float型別。如果使用這個型別，儲存的資料範圍是沒有什麼大的問題的，我們在網路上看到的文章大部分也是使用這個型別來儲存結果的。但是，我在實踐中多次遇到用float型別得不到正確的結果的問題，後來發現核心的原因是float的計算精度嚴重不足，特別是對於積分圖這種連續的加法的計算，累計誤差會越來越嚴重，當計算量足夠大時，就會出現明顯的誤差瑕疵。因此，這個資料型別從本質上來說，對積分圖是不夠安全的。

　　關於這一點，實際上已經有不少作者注意到了，我在博文：SSE影像演算法最佳化系列十四：區域性均方差及區域性平方差演算法的最佳化中也有提及。

　　第三個是double型別，這個型別也是64位的，因此，資料範圍毫無問題，計算精度經過測試，也是沒有什麼問題的，而且，編譯器和指令集的支援和最佳化做的都還很不錯，因此，個人認為這個資料型別是用來儲存積分圖最為合適的型別，但是有一個不友好的特點，計算速度慢，而且指令集對其能加速的空間有限。

　　2、積分圖雖然能做到某些演算法和引數無關，但是其並不是最佳的最速度

　　使用積分圖技術，首先是要分配積分圖佔用的那部分額外的記憶體（而且是相當客觀的記憶體），其次，積分圖本身的計算也是需要一定時間的。還有，積分圖的計算必須對邊界部分做特別的判斷，這個當引數較大時，計算量有的時候還是相當可觀的。

　　還是回到我們的非區域性均值濾波上吧。上面說了這麼多，意思就是雖然非區域性均值可以用積分圖去最佳化，但是還是不是很好，那有沒有更好的更快的實現呢，其實在我的下面兩篇部落格裡就已經有了相關的技術。

　　一是： SSE影像演算法最佳化系列十三：超高速BoxBlur演算法的實現和最佳化（Opencv的速度的五倍）

二是： SSE影像演算法最佳化系列十四：區域性均方差及區域性平方差演算法的最佳化

　　因為非區域性均值裡使用的積分圖是差的平方的積分圖，因此我們只要對上述參考博文一里面的參與積分的資料稍微換一下即可，一個簡單的C++程式碼如下所示：

//    SearchRadius    搜尋的領域半徑
//    PatchRadius      計算相似度的塊的半徑
//    Delta     高斯平滑的引數
int IM_NLM_Denoising(unsigned char *Src, unsigned char *Dest, int Width, int Height, int Stride, int SearchRadius, int PatchRadius, float Delta)
{
    int Status = IM_STATUS_OK;
    int Channel = Stride / Width;
    if (Src == NULL)                        return IM_STATUS_NULLREFRENCE;
    if ((Width <= 0) || (Height <= 0))        return IM_STATUS_INVALIDPARAMETER;
    if ((Channel != 1))                        return IM_STATUS_INVALIDPARAMETER;    

    int ExpandW = Width + SearchRadius + PatchRadius + SearchRadius + PatchRadius;
    int ExpandH = Height + SearchRadius + PatchRadius + SearchRadius + PatchRadius;
    int DiffW = Width + PatchRadius + PatchRadius;
    int DiffH = Height + PatchRadius + PatchRadius;

    float *Weight = (float *)calloc(Width * Height , sizeof(float));
    float *Sum = (float *)calloc(Width * Height , sizeof(float));
    int *ColValue = (int *)malloc((Width + PatchRadius + PatchRadius) * sizeof(int));
    unsigned char *Expand = (unsigned char *)malloc(ExpandH * ExpandW * sizeof(unsigned char));
    if ((Weight == NULL) || (Sum == NULL) || (ColValue == NULL) || (Expand == NULL))
    {
        Status = IM_STATUS_OUTOFMEMORY;
        goto FreeMemory;
    }
    Status = IM_GetExpandImage(Src, Expand, Width, Height, Stride, ExpandW, SearchRadius + PatchRadius, SearchRadius + PatchRadius, SearchRadius + PatchRadius, SearchRadius + PatchRadius, IM_EDGE_MIRROR);
    if (Status != IM_STATUS_OK)    goto FreeMemory;

    int Area = (2 * PatchRadius + 1) * (2 * PatchRadius + 1);
    float Inv = -1.0f / Area / Delta / Delta;


    for (int YY = -SearchRadius; YY <= SearchRadius; YY++)
    {
        for (int XX = -SearchRadius; XX <= SearchRadius; XX++)
        {
            for (int Y = 0; Y < Height; Y++)
            {
                float *LinePS = Sum + Y * Width;
                float *LinePW = Weight + Y * Width;
                unsigned char *LinePE = Expand + (Y + YY + SearchRadius + PatchRadius) * ExpandW + XX + SearchRadius + PatchRadius;
                if (Y == 0)
                {
                    memset(ColValue, 0, DiffW * sizeof(int));
                    for (int Z = -PatchRadius; Z <= PatchRadius; Z++)
                    {
                        unsigned char *LineP1 = Expand + (Z + PatchRadius + YY + SearchRadius) * ExpandW + XX + SearchRadius;
                        unsigned char *LineP2 = Expand + (Z + PatchRadius + SearchRadius) * ExpandW + SearchRadius;
                        for (int X = 0; X < DiffW; X++)
                        {
                            int Value = LineP2[X] - LineP1[X];
                            ColValue[X] += Value * Value;
                        }
                    }
                }
                else
                {
                    unsigned char *LineOut1 = Expand + (Y - 1 + YY + SearchRadius) * ExpandW + XX + SearchRadius;
                    unsigned char *LineOut2 = Expand + (Y - 1 + SearchRadius) * ExpandW + SearchRadius;
                    unsigned char *LineIn1 = Expand + (Y + PatchRadius + PatchRadius + YY + SearchRadius) * ExpandW + XX + SearchRadius;
                    unsigned char *LineIn2 = Expand + (Y + PatchRadius + PatchRadius + SearchRadius) * ExpandW + SearchRadius;
                    for (int X = 0; X < DiffW; X++)
                    {
                        int Out = LineOut2[X] - LineOut1[X];
                        int In = LineIn2[X] - LineIn1[X];
                        ColValue[X] -= Out * Out - In * In;                                    //    更新列資料
                    }
                }
                int SumA = IM_SumofArray(ColValue, PatchRadius * 2 + 1);                //    處理每行第一個資料    
                float W = IM_Fexp(SumA * Inv);
                LinePW[0] += W;
                LinePS[0] += W * LinePE[0];
                for (int X = 1; X < Width; X++)
                {
                    SumA = SumA - ColValue[X - 1] + ColValue[X + PatchRadius + PatchRadius];
                    float W = IM_Fexp(SumA * Inv);
                    LinePW[X] += W;
                    LinePS[X] += W * LinePE[X];
                }
            }
        }
    }
    for (int Y = 0; Y < Height; Y++)
    {
        int Index = Y * Width;
        unsigned char *LinePD = Dest + Y * Stride;
        for (int X = 0; X < Width; X++)
        {
            LinePD[X] = Sum[Index + X] / (Weight[Index + X]);
        }
    }
FreeMemory:
    if (Expand != NULL)            free(Expand);
    if (Weight != NULL)            free(Weight);
    if (Sum != NULL)               free(Sum);
    if (ColValue != NULL)         free(ColValue);
    return Status;
}

　　對於常用的搜尋半徑為10，塊半徑為3, 500 * 500的灰度圖耗時大概是500ms，相當的慢啊。

　　透過我一系列博文裡的資料，可以知道上面的迴圈內的程式碼其實是很容易進行指令集最佳化的，基本上我那個方框模糊的最佳化是同一個技巧，經過指令集最佳化後，500*500的灰度圖的耗時大概在200ms左右，如果加上執行緒技術，可以最佳化到75ms左右，這個時間和非區域性均值濾波（NL-means）演算法的CUDA最佳化加速文章裡提的CUDA的最佳化速度基本已經差不多了。

　　另外，在搜尋半徑較小時，一種可行的最佳化方式是進行行列分離的卷積，即先計算中心行的結果，然後已這個結果為原始資料，在計算列方向的卷積，注意不是同時計算中心列和中心行的累加，而是有順序的，這樣就可以利用到周邊領域所有的資訊，這個時候計算量就會大為下降，計算的耗時也可以明顯提高。這種最佳化方式必須測試是否對去燥的效果有很大的影響，因為他有可能會產生較為明顯的水平或垂直線條效果。

　　另外，如果搜尋半徑較大，還可以嘗試在上述基礎上再進行45度和135度兩個方向的卷積，以便抵消這種線條效果，同樣提速也還是很明顯的。

　　作為一個特例，有些情況下我們可能需要搜尋半徑為2，塊大小也為2的非區域性均值濾波，這種尺寸的濾波對於高強度的高斯噪音是沒有什麼去燥效果的，但是對於小規模的噪音，比如一般視訊會議中的噪音還是有較強的抑制作用，因此，也還是有應用場景的，但是這種場景對演算法的速度提出了極高的要求，如果考慮流暢性，給這個演算法的處理不易超過20ms，而現在的影片流越來越高畫質了，因此，對這個演算法的最佳化處理就必須做的更到位。

　　針對這個特例，我又做了一些最佳化，首先，因為是小半徑，所以可以使用行列分離的演算法，即先計算如下區域的合成結果：

　　計算得到中間的結果後再透過中間結果計算下述區域的合成結果作為最終值：

　　這裡面的最佳化技巧有很多很多。

　　我這裡提幾個想法供參考：　　

　　1、中心點C的處的權重不需要計算，因為必然為1。

　　2、在進行水平計算時，去除掉中線後，只有4個點了，PP/P/N/NN，我們可以透過一些組合手段把他們的權重以及權重乘以值的量一次性計算出來，這樣就可以直接把中間值給搞定了，如此做的好處是，不用多次儲存累加的權重值和乘以值，也就是說捨棄掉了前面程式碼裡的Weight和Sum變數，這對演算法速度的提高也是有極大的貢獻的，垂直的計算也是一樣的。

3、那個exp函式是個比較耗時的函式，在http://martin.ankerl.com/2007/10/04/optimized-pow-approximation-for-java-and-c-c/中一個比較快速的函式，雖然精度一般，但是經過測試，可以滿足本演算法的需求，但是其計算量小很多。

　　4、去掉中心點後，無論在垂直或者水平方向都只有4個資料，這為SIMD指令的某些讀取和最佳化提供了無尚的便利和方便，而有些資料真的是巧合加天工設計，比如水平方向處理時，我們需要一次性處理4個點（出去中間那個），而我們用_mm_loadl_epi64恰好可以讀取8個位元組，這個時候，讀取的8個位元組如下：

　　　　　　　　A0 　　A1 　　A2 　　A3 　　A4 　　A5 　　A6 　　A7

　　正好可以組成4對5個組合 :

　　　　　　　　　　A0　　A1　　A2　　A3　　A4

　　　　　　　　　　A1 　　A2 　 A3 　　A4 　　A5

　　　　　　　　　　A2 　　A3 　A4 　　A5 　　A6

　　　　　　　　　　A3 　　A4 　A5 　　A6 A7

　　如果多一個位元組，都不好處理了（主要是不好讀取），多麼完美的事情啊。

　　透過多種最佳化方式後，對於常見的1080P影片（1920*1080），處理其Y分量（灰度）耗時能做到28ms了，如果使用2個執行緒或4個執行緒，則可以完全滿足實時的需求，如果是720P的影片，則單核也能到20ms。