導向濾波算法的實現

由上篇導向濾波算法分析，根據(5)~(8)式就能夠計算輸出圖像Q

　　(5)

　　(6)

　　(7)

　　(8)

其中，，/a_i和/b_i的結果要計算全部覆蓋了像素i的窗口W_k的a_k和b_k的平均值。除了用平均值，在實際應用中，我還看到過其餘的計算/a_i和/b_i的方法。好比根據像素i在窗口W_k的位置，給予不一樣的權重。若是i距離窗口W_k的中心位置越遠，則給予a_k和b_k越低的權重。若是i位於窗口中心，則a_k和b_k有最高的權重。最經常使用的就是用高斯分佈來給予不一樣的權重。這樣考慮到了附近像素距離遠近對i影響的大小，最後的結果會更精確一些。

這裏我仍是用最簡單的平均值的方法來計算/a_i和/b_i。咱們前面已經假定了在窗口W_k內，a_k和b_k是常數，所以a_k和b_k只和W_k的位置有關。取W_k爲半徑爲r的方形窗口，使窗口的中心像素位置遍歷整個圖像，那麼就使W_k取到了不一樣的全部位置，在每一個位置計算出相應的a_k和b_k。全部的a_k和b_k組成了和輸入圖像P相同長寬維度的數據集合，記爲A和B。對於任意像素i，/a_i和/b_i即分別爲以i爲中心半徑r的窗口W_k內A和B的數據均值，這不正是咱們熟悉的圖像均值模糊的算法嗎？而圖像均值模糊有幾種比較成熟的快速算法，好比積分圖算法和基於直方圖的快速算法。只要有了A和B，就能夠方便的應用均值模糊得出/a_i和/b_i，從而應用(8)計算出輸出圖像Q。

爲了計算A和B，從(6)式看到，須要分別計算輸入圖像P和導向圖G的均值模糊結果。而(5)式須要計算導向圖G的方差，還有P和G的協方差。方差和協方差都涉及到乘積的求和計算，能夠由下面的公式，經過積分圖來計算。這兩個公式很容易推導出來，就不贅述了。

一個平方和，一個乘積和均可以用積分圖來計算。只是要注意當圖像足夠大的時候，要用合適的數據類型。假設像素的數據範圍是0~255的整型，若是平方積分圖用32位的整型數據，那麼只能支持最大256x256大小的圖像。超過這個大小，就必需要用64位的整型了。下面給出使用模板函數的乘積積分圖函數，能夠根據須要使用不一樣的數據類型。p1和p2是圖像數據指針，當它們指向相同的數據時，這個函數就變成了平方積分圖。注意積分圖的長寬比原始數據的長寬都要大1。

/* Cumulative image of the multiplication of p1.*p2. p1 and p2 can be the same pointer and it becomes square cumulative image. The returned cumulative image MUST be freed by the caller! */ template <class T1, class T2, class T3> BOOL MultiplyCumImage(T1 *p1, T2 *p2, T3 **ppCum) { long i, j; long Width = GetWidth(); long Height = GetHeight(); long integral_len = (Width + 1)*(Height + 1); // Only allocate cumulative image memory when *ppCum is NULL
    if (*ppCum == NULL) { try { *ppCum = new T3[integral_len]; } catch (CException *pe) { pe->Delete(); return FALSE; } } memset(*ppCum, 0, sizeof(T3)*integral_len); // The cumulative values of the leftmost and the topmost pixels are always 0.
    for (i = 1; i <= Height; i++) { T3 *prow, *puprow; prow = *ppCum + i*(Width + 1); puprow = *ppCum + (i - 1)*(Width + 1); T3 sum = 0; long up_row_idx = (i - 1)*Width; for (j = 1; j <= Width; j++) { long idx = up_row_idx + j - 1; sum += p1[idx] * p2[idx]; prow[j] = puprow[j] + sum; } } return TRUE; }

View Code

 這樣導向濾波實現的主要問題都解決了，算法步驟以下：

1. 計算引導圖G的積分圖和平方積分圖
2. 計算輸入圖像P的積分圖， P和G的乘積積分圖
3. 用上兩步得出的積分圖計算P和G的均值，G的方差，P和G的協方差，窗口半徑爲r
4. 而後用(5)(6)式計算係數圖A和B
5. 計算A和B的積分圖
6. 計算A和B的窗口半徑r的均值，並用(8)式計算輸出圖Q

 下面的參考代碼中，pData存儲輸入和輸出圖像，pGuidedData引導圖，radius領域半徑

    long len = Width*Height; // Cululative image and square cululative for guided image G
    UINT64 *pCum = NULL, *pSquareCum = NULL; CumImage(pGuidedData, &pCum); MultiplyCumImage(pGuidedData, pGuidedData, &pSquareCum); // Allocate memory for a and b
    float *pa, *pb; pa = new float[len]; pb = new float[len]; memset(pa, 0, sizeof(float)*len); memset(pb, 0, sizeof(float)*len); UINT64 *pInputCum = NULL, *pGPCum = NULL; CumImage(pData, &pInputCum); MultiplyCumImage(pGuidedData, pData, &pGPCum); int field_size; UINT64 cum, square_cum; long uprow, downrow, upidx, downidx;        // In cumulative image
    long leftcol, rightcol; float g_mean, g_var;        // mean and variance of guided image
    long row_idx = 0; UINT64 p_cum, gp_cum; float p_mean; // Calculate a and b // Since we're going to calculate cumulative image of a and b, we have to calculate the whole image of a and b.
    for (i = 0; i < Height; i++) { // Check the boundary for radius
        if (i < radius) uprow = 0; else uprow = i - radius; upidx = uprow*(Width + 1); if (i + radius >= Height) downrow = Height; else downrow = i + radius + 1; downidx = downrow*(Width + 1); for (j = 0; j < Width; j++) { // Check the boundary for radius
            if (j < radius) leftcol = 0; else leftcol = j - radius; if (j + radius >= Width) rightcol = Width; else rightcol = j + radius + 1; field_size = (downrow - uprow)*(rightcol - leftcol); long p1, p2, p3, p4; p1 = downidx + rightcol; p2 = downidx + leftcol; p3 = upidx + rightcol; p4 = upidx + leftcol; // Guided image summary in the field
            cum = pCum[p1] - pCum[p2] - pCum[p3] + pCum[p4]; // Guided image square summary in the field
            square_cum = pSquareCum[p1] - pSquareCum[p2] - pSquareCum[p3] + pSquareCum[p4]; // Field mean
            g_mean = (float)(cum) / field_size; // Field variance
            g_var = float(square_cum) / field_size - g_mean * g_mean; // Summary of input image in the field
            p_cum = pInputCum[p1] - pInputCum[p2] - pInputCum[p3] + pInputCum[p4]; // Input image field mean
            p_mean = float(p_cum) / field_size; // Multiply summary in the field
            gp_cum = pGPCum[p1] - pGPCum[p2] - pGPCum[p3] + pGPCum[p4]; long idx = row_idx + j; pa[idx] = (float(gp_cum) / field_size - g_mean*p_mean) / (g_var + epsilon); pb[idx] = p_mean - g_mean*pa[idx]; } row_idx += Width; } // not needed after this
    delete[] pCum; delete[] pSquareCum; delete[] pInputCum; delete[] pGPCum; // Cumulative image of a and b
    float *pCuma = NULL, *pCumb = NULL; CumImage(pa, &pCuma); CumImage(pb, &pCumb); // Finally calculate the output image q=ag+b
    float mean_a, mean_b; row_idx = Hstart*Width; for (i = Hstart; i < Hend; i++) { // Check the boundary for radius
        if (i < radius) uprow = 0; else uprow = i - radius; upidx = uprow*(Width + 1); if (i + radius >= Height) downrow = Height; else downrow = i + radius + 1; downidx = downrow*(Width + 1); for (j = Wstart; j < Wend; j++) { // Check the boundary for radius
            if (j < radius) leftcol = 0; else leftcol = j - radius; if (j + radius >= Width) rightcol = Width; else rightcol = j + radius + 1; field_size = (downrow - uprow)*(rightcol - leftcol); long p1, p2, p3, p4; p1 = downidx + rightcol; p2 = downidx + leftcol; p3 = upidx + rightcol; p4 = upidx + leftcol; // Field mean
            mean_a = (pCuma[p1] - pCuma[p2] - pCuma[p3] + pCuma[p4]) / field_size; // Field mean
            mean_b = (pCumb[p1] - pCumb[p2] - pCumb[p3] + pCumb[p4]) / field_size; // New pixel value
            long idx = row_idx + j; int value = int(mean_a*pGuidedData[idx] + mean_b); CLAMP0255(value); pData[idx] = value; } row_idx += Width; } delete[] pa; delete[] pb; delete[] pCuma; delete[] pCumb;

View Code

 導向濾波還有一種快速算法，基本思想是經過下采樣輸入圖P和引導圖G，獲得較小的圖像P'和G'。用它們來計算係數A'和B'。而後經過線性插值的辦法恢復原始大小獲得近似的A和B，用來和原始大小的引導圖來計算輸出Q。這樣，a_k和b_k的計算不是在全尺寸圖像上，能節省很多運算量，而最終的結果不受很大影響。

 下面以帶有保邊平滑特性的導向濾波爲例，來看看效果。如上一篇所說，當輸入圖P和引導圖G相同時，導向濾波呈現保邊平滑特性。

  

                              原圖                                                                                                            半徑3，ε=5²

 

                             半徑3，ε=10²                                                                         半徑3，ε=20²

  

                原圖                                                         半徑5，ε=20²

 整體來講，P=G時，導向濾波的保邊平滑特性和帶有保邊功能領域平滑濾波有相似的效果。