O(1)效率的表面模糊算法優化。

通過各種手段的優化，實作O(1)效率的表面模糊效果，并共享全部C++源代碼。

很久沒有寫文章了，主要是最近一段時間沒有以前那麼多空暇空間，記憶體和CPU占用率一緻都很高，應前幾日群裡網友的要求，今天發個表面模糊的小程式來找回之前寫部落格的熱情吧。

國内我認為，破解表面模糊的原理的最早作者是我一直很崇拜的一位女士，她不會程式設計，英文也不怎麼好，僅憑電腦和Excel兩個工具破解了PS了很多算法，真是個巾帼英雄。

詳見位址：http://www.missyuan.com/thread-428384-1-1.htm

網上的有關該算法的matlab實作參考：http://www.cnblogs.com/tiandsp/archive/2012/11/06/2756441.html

用C實作的參考：http://blog.csdn.net/maozefa/article/details/8270990

表面模糊是屬于典型的EPF濾波器中的一種，在PS的架構下好像也隻有這一種自帶的EPF算法，其核心也是數卷積的範疇，隻是卷積的核是随着内容而變的，也屬于方形半徑内的算法，借助于直方圖是可以做到于參數無關的O(1)算法。關于直方圖的相關架構參考我的博文：任意半徑局部直方圖類算法在PC中快速實作的架構。，但本文代碼對其做了稍許改動。

為了表述友善，我們以灰階圖像為例進行說明。首先，表面模糊有兩個參數，半徑Radius和門檻值Threshold。如果我們知道了以某點為中心，半徑為Radius範圍内的直方圖資料Hist，以及該點的像素值，那根據原始的算法，其計算公式為：

//　　最原始的算法
void Calc(unsigned short *Hist, unsigned char Value, int Threshold, unsigned char *&Pixel)
{
    int Weight, Sum = 0, Divisor = 0;
    for (int Y = 0; Y < 256; Y++)
    {
        Weight = Hist[Y] * (2500 - abs(Y - Value) * 1000 / Threshold);
        if (Weight < 0) Weight = 0;
        Sum += Weight * Y;
        Divisor += Weight;
    }
    if (Divisor > 0) *Pixel = (Sum + (Divisor >> 1)) / Divisor;
}

　注意這裡我們為了減少浮點計算，将權重的計算公式放大了2500倍以便進行定點化，同時必須在最後增加一個Divisor > 0的判斷，因為當Threshold很小時，可能會出現Divisor為0的現象。

上述代碼針對1000*1000的灰階圖的執行時間約為1250ms,其中直方圖的更新時間隻有約50ms，速度難以接受。

分析計算方法1，很明顯權重計算的幾個加減乘除以及下面的那句判斷是比較耗時的，而其隻是Y-Value的一個函數，是以，我們可以提前建立一個表，該表的索引範圍從Min[Y - Value]到Max[Y - Value]之間，很明顯，這個範圍是[-255, 255]，是以，建立如下的一個查找表：

for (int Y = -255; Y <= 255; Y++)
{
    int Factor = (2500 - abs(Y) * 1000 / Threshold);       
    if (Factor < 0) Factor = 0;
    Intensity[Y + 255] = Factor;
}

　　有了這個查找表，我們來實作第二個版本的算法如下：

//    改進後的算法
unsigned char Calc2(unsigned short *Hist, unsigned char Value, unsigned short *Intensity)
{
    int Weight = 0, Sum = 0, Divisor = 0;
    unsigned short *Offset = Intensity + 255 - Value;
    for (int Y = 0; Y < 256; Y++)
    {
        Weight = Hist[Y] * Offset[Y];
        Sum += Weight * Y;
        Divisor += Weight;
    }
    if (Divisor > 0)
        return (Sum + (Divisor >> 1)) / Divisor;        //    四舍五入
    else
        return Value;
}

　　同樣大小的圖，執行時間為350ms，速度提高約為3倍。

我們接着來思考問題，上述有256個循環，如果我們将循環手動展開，會不會有提高呢，我們先把代碼更改如下：

//    優化後的算法
unsigned char Calc3(unsigned short *Hist, unsigned char Value, unsigned short *Intensity)       
{
    int Weight = 0, Sum = 0, Divisor = 0;
    unsigned short *Offset = Intensity + 255 - Value;
    Weight = Hist[0] * Offset[0];
    Sum += Weight * 0;  Divisor += Weight;        //    能不能用使用指令集的并行，沒有去測試了
    Weight = Hist[1] * Offset[1];
    Sum += Weight * 1;  Divisor += Weight;
    Weight = Hist[2] * Offset[2];
    Sum += Weight * 2;  Divisor += Weight;
    Weight = Hist[3] * Offset[3];
    Sum += Weight * 3;  Divisor += Weight;
   /////////////////////////// ............................................................................
    Weight = Hist[251] * Offset[251];
    Sum += Weight * 251;  Divisor += Weight;
    Weight = Hist[252] * Offset[252];
    Sum += Weight * 252;  Divisor += Weight;
    Weight = Hist[253] * Offset[253];
    Sum += Weight * 253;  Divisor += Weight;
    Weight = Hist[254] * Offset[254];
    Sum += Weight * 254;  Divisor += Weight;
    Weight = Hist[255] * Offset[255];
    Sum += Weight * 255;  Divisor += Weight;
    if (Divisor > 0)
        return (Sum + (Divisor >> 1)) / Divisor;        //    四舍五入
    else
        return Value;
}

　　為表述友善，中間省略了一些代碼。

測試結果為250ms，又快了一點點，為什麼呢，我分析認為第一是減少了循環計數的時間，第二循環展開的乘以常數會被CPU優化為相關的移位或其他操作，而Calc2内部編譯器是無法優化的。

這樣的函數系統一般是不會内聯的，即使你在函數前面加上inline辨別符，但是你可以在前面加上__forceinline辨別，強制他内聯，但是如果你這樣做，你會發現速度反而會嚴重下降，為什麼，請大家自行分析。

我們在自己仔細看看，上面的循環很容易用SSE函數實作，既然我們的直方圖的擷取和更新利用了SSE，這裡為什麼不用呢，這樣就誕生了我們的Calc4函數。

//    用SSE優化的算法
unsigned char Calc4(unsigned short *Hist, unsigned char Value, unsigned short *Intensity, unsigned short *Level)
{
    unsigned short *Offset = Intensity + 255 - Value;
    __m128i SumS = _mm_setzero_si128();
    __m128i WeightS = _mm_setzero_si128();
    for (int K = 0; K < 256; K += 8)
    {
        __m128i H = _mm_load_si128((__m128i const *)(Hist + K));
        __m128i L = _mm_load_si128((__m128i const *)(Level + K));                //    有能力可以使用256位的AVX寄存器
        __m128i I = _mm_loadu_si128((__m128i const *)(Offset + K));
        SumS = _mm_add_epi32(_mm_madd_epi16(_mm_mullo_epi16(L, I), H), SumS);
        WeightS = _mm_add_epi32(_mm_madd_epi16(H, I), WeightS);
    }
    const int *WW = (const int *)&WeightS;
    const int *SS = (const int *)&SumS;

    int Sum = SS[0] + SS[1] + SS[2] + SS[3];
    int Divisor = WW[0] + WW[1] + WW[2] + WW[3];
    if (Divisor > 0)
        return (Sum + (Divisor >> 1)) / Divisor;        //    四舍五入
    else
        return Value;
}

　　關于上面幾個SSE函數的使用，我不想多談，也沒啥難易了解的，注意其中的Level是我們為了友善，預定義的一個表，其形式如下：

for (int Y = 0; Y < 256; Y++)    Level[Y] = Y;            //    這個是為CalcSSE友善的使用的，其他兩可以删除掉這裡

不定義這個也應該可以由其他的SSE函數構造k/k+1/k+2/k+3/k+4/k+5/k+6/k+7這樣的__m128i變量，我這裡這樣做隻是為了友善，你也可以自己更改下。

我們直接把Calc4嵌入到程式中，運作，發現運作時間降低到了100ms，比Calc3有提高了2倍多，但是效果似乎不對，怎麼回事呢。

這主要是因為上述的SSE函數是針對unsigned short類型，而我們構造的Intensity資料較大，進行乘法後會超出unsigned short所能表達的範圍，是以我們需要改動Intensity的定義：

//    為了SSE裡不溢出，把這裡的資料變小，當然這樣算法的準确度降低了，但是為了速度.......
    for (int Y = -255; Y <= 255; Y++)
    {
        int Factor = (255 - abs(Y) * 100 / Threshold);        
        if (Factor < 0) Factor = 0;
        Intensity[Y + 255] = Factor / 2;
    }

　　最後一個除以2估計是因為SSE内部還是按照signed short處理的，這樣做會導緻算法的精度降低。

經過上述改動，效果就正确了。

對于彩色圖像，一種做法就是直接擴充現在單通道的代碼，讓其支援三通道，另外一個辦法就是把圖像先拆分成3通道獨立的資料，然後沒通道獨立處理，處理完成後再合成，這樣做有兩個好處，第一是代碼複用；第二就是如果支援Openmp或者其他的并行庫，可以讓3通道并行起來執行。但是也有2個不足，第一是記憶體占用會增加很多，因為這種算法是不支援In-Place操作的，是以必須配置設定6份單通道的資料，而算法内部配置設定的記憶體由于并行的關系也要增加一些（不是三倍），及時考慮到可以把其中三個通道的數放置到Dest中，也會增加3份通道的資料，這對于某些裝置可能是難以接受的（比如低端的安卓機)。具體如何使用就看應用場景了。

針對實際的應用，一種可選的進一步加速的方式就是把圖像的色階範圍進一步縮小，比如由256色階變為128或者64色階，這樣理論上還可以在快2倍到4倍，不過效果會稍有下降，一般128位時還是可以接受的。

本文的完整VS2013代碼下載下傳位址（解壓密碼本人部落格名）：https://files.cnblogs.com/files/Imageshop/SurfaceBlur.rar

我看到很多人轉載我的文章，我很感謝，但是很多人沒有一點點的尊重别人的意識，轉載請你在博文的最前面聲明為轉載，并不要更改本文下部打賞二維碼。

****************************作者： laviewpbt 時間： 2015.10.24 聯系QQ: 33184777 轉載請保留本行資訊**********************

O(1)效率的表面模糊算法優化。

繼續閱讀

C++圖像處理 -- 表面模糊