綜述：圖像濾波常用算法實現(xiàn)及原理解析（2）

雙邊濾波的代碼

opencv中提供了bilateralFilter()函數(shù)來實現(xiàn)雙邊濾波操作，其原型如下：

void cv::bilateralFilter(InputArray src,
OutputArray  dst,
int  d,
double  sigmaColor,
double  sigmaSpace,
int  borderType = BORDER_DEFAULT 
) 

• InputArray src: 輸入圖像，可以是Mat類型，圖像必須是8位整型或浮點型單通道、三通道的圖像。
• OutputArray dst: 輸出圖像，和原圖像有相同的尺寸和類型。
• int d: 表示在過濾過程中每個像素鄰域的直徑范圍。如果這個值是非正數(shù)，則函數(shù)會從第五個參數(shù)sigmaSpace計算該值。
• double sigmaColor: 顏色空間過濾器的值，這個參數(shù)的值月大，表明該像素鄰域內(nèi)有越寬廣的顏色會被混合到一起，產(chǎn)生較大的半相等顏色區(qū)域。 （這個參數(shù)可以理解為值域核的 和 )
• double sigmaSpace: 坐標空間中濾波器的sigma值，如果該值較大，則意味著越遠的像素將相互影響，從而使更大的區(qū)域中足夠相似的顏色獲取相同的顏色。當d>0時，d指定了鄰域大小且與sigmaSpace無關(guān)，否則d正比于sigmaSpace. （這個參數(shù)可以理解為空間域核的 和 ）
• int borderType=BORDER_DEFAULT: 用于推斷圖像外部像素的某種邊界模式，有默認值BORDER_DEFAULT.

具體代碼如下：

#include <iostream>#include <opencv2/opencv.hpp>
using namespace std;using namespace cv;
//定義全局變量const int g_ndMaxValue = 100;const int g_nsigmaColorMaxValue = 200;const int g_nsigmaSpaceMaxValue = 200;int g_ndValue;int g_nsigmaColorValue;int g_nsigmaSpaceValue;
Mat g_srcImage;Mat g_dstImage;
//定義回調(diào)函數(shù)void on_bilateralFilterTrackbar(int, void*);
int main(){    g_srcImage = imread("lena.jpg");
    //判斷圖像是否加載成功    if(g_srcImage.empty())    {        cout << "圖像加載失敗!" << endl;        return -1;    }    else        cout << "圖像加載成功!" << endl << endl;
    namedWindow("原圖像", WINDOW_AUTOSIZE);    imshow("原圖像", g_srcImage);
    //定義輸出圖像窗口屬性和軌跡條屬性    namedWindow("雙邊濾波圖像", WINDOW_AUTOSIZE);    g_ndValue = 10;    g_nsigmaColorValue = 10;    g_nsigmaSpaceValue = 10;
    char dName[20];    sprintf(dName, "鄰域直徑 %d", g_ndMaxValue);
    char sigmaColorName[20];    sprintf(sigmaColorName, "sigmaColor %d", g_nsigmaColorMaxValue);
    char sigmaSpaceName[20];    sprintf(sigmaSpaceName, "sigmaSpace %d", g_nsigmaSpaceMaxValue);
    //創(chuàng)建軌跡條    createTrackbar(dName, "雙邊濾波圖像", &g_ndValue, g_ndMaxValue, on_bilateralFilterTrackbar);    on_bilateralFilterTrackbar(g_ndValue, 0);
    createTrackbar(sigmaColorName, "雙邊濾波圖像", &g_nsigmaColorValue,                     g_nsigmaColorMaxValue, on_bilateralFilterTrackbar);    on_bilateralFilterTrackbar(g_nsigmaColorValue, 0);
    createTrackbar(sigmaSpaceName, "雙邊濾波圖像", &g_nsigmaSpaceValue,                    g_nsigmaSpaceMaxValue, on_bilateralFilterTrackbar);    on_bilateralFilterTrackbar(g_nsigmaSpaceValue, 0);
    waitKey(0);
    return 0;}
void on_bilateralFilterTrackbar(int, void*){    bilateralFilter(g_srcImage, g_dstImage, g_ndValue, g_nsigmaColorValue, g_nsigmaSpaceValue);    imshow("雙邊濾波圖像", g_dstImage);}

需要有高斯濾波和雙邊濾波的相關(guān)知識背景才能更好的理解導(dǎo)向濾波。在導(dǎo)向濾波中，首先利用了局部線性模型。這個模型認為某函數(shù)上一點與其近鄰部分的點成線性關(guān)系，一個復(fù)雜的函數(shù)就可以用很多局部的線性函數(shù)來表示，當需要求該函數(shù)上某一點的值時，只需要計算所有包含該點的線性函數(shù)的值并取平均值即可。這種模型，在表示非解析函數(shù)上，非常有用。

同理，我們可以認為圖像是一個二維函數(shù)，并且假設(shè)該函數(shù)的輸出與輸入在一個二維窗口內(nèi)滿足線性關(guān)系，如下：

其中，是輸出像素的值，是輸入圖像的值，和是像素索引，和是當窗口中心位于k時該線性函數(shù)的系數(shù)。其實，輸入圖像不一定是待濾波的圖像本身，也可以是其他圖像即引導(dǎo)圖像，這也是為何稱為引導(dǎo)濾波的原因。對上式兩邊取梯度，可以得到:

即當輸入圖像有梯度時，輸出也有類似的梯度，現(xiàn)在可以解釋為什么引導(dǎo)濾波有邊緣保持特性了。下一步是求出線性函數(shù)的系數(shù)，也就是線性回歸，即希望擬合函數(shù)的輸出值與真實值之間的差距最小，也就是讓下式最小:

這里只能是待濾波圖像，并不像那樣可以是其他圖像。同時，之前的系數(shù)用于防止求得的過大，也是調(diào)節(jié)濾波器濾波效果的重要參數(shù)(相當于L2正則化的權(quán)重懲罰)。接下來利用最小二乘法的原理令  和 得到2個二元一次方程，求解得到：

其中 是 在窗口的平均值， 是 在窗口 的方差， 是窗口 中的像素個數(shù)，  是待濾波圖像在窗口 中的均值。在計算每個窗口的線性系數(shù)時，我們可以發(fā)現(xiàn)一個像素會被多個窗口包含，也就是說，每個像素都由多個線性函數(shù)所描述。因此，如之前所說，要具體求某一點的輸出值時，只需將所有包含該點的線性函數(shù)值平均即可，如下:

這里， 是所有包含像素 的窗口， 是其中心位置。

當把引導(dǎo)濾波用作邊緣保持濾波器時，往往有  ，如果 ，顯然是為最小值的解，從上式可以看出，這時的濾波器沒有任何作用，將輸入原封不動的輸出。如果 ，在像素強度變化小的區(qū)域（或單色區(qū)域），有近似于（或等于0，而近似于（或等于） ，即做了一個加權(quán)均值濾波；而在變化大的區(qū)域，近似于1，近似于0，對圖像的濾波效果很弱，有助于保持邊緣。而 的作用就是界定什么是變化大，什么是變化小。在窗口大小不變的情況下，隨著的增大，濾波效果越明顯。

在濾波效果上，引導(dǎo)濾波和雙邊濾波差不多，然后在一些細節(jié)上，引導(dǎo)濾波較好(在PS的磨皮美白中，經(jīng)過實踐，效果更好)。引導(dǎo)濾波最大的優(yōu)勢在于，可以寫出時間復(fù)雜度與窗口大小無關(guān)的算法，因此在使用大窗口處理圖片時，其效率更高。

同樣，OpenCV中也有導(dǎo)向濾波的接口。具體代碼如下：

void cv::ximgproc::guidedFilter ( InputArray  guide,
InputArray  src,
OutputArray  dst,
int  radius,
double  eps,
int  dDepth = -1 
) 

這邊有個基于scipy實現(xiàn)的python代碼，可以參考一下：

import numpy as npimport scipy as spimport scipy.ndimage

def box(img, r):    """ O(1) box filter        img - >= 2d image        r   - radius of box filter    """    (rows, cols) = img.shape[:2]    imDst = np.zeros_like(img)

    tile = [1] * img.ndim    tile[0] = r    imCum = np.cumsum(img, 0)    imDst[0:r+1, :, ...] = imCum[r:2*r+1, :, ...]    imDst[r+1:rows-r, :, ...] = imCum[2*r+1:rows, :, ...] - imCum[0:rows-2*r-1, :, ...]    imDst[rows-r:rows, :, ...] = np.tile(imCum[rows-1:rows, :, ...], tile) - imCum[rows-2*r-1:rows-r-1, :, ...]
    tile = [1] * img.ndim    tile[1] = r    imCum = np.cumsum(imDst, 1)    imDst[:, 0:r+1, ...] = imCum[:, r:2*r+1, ...]    imDst[:, r+1:cols-r, ...] = imCum[:, 2*r+1 : cols, ...] - imCum[:, 0 : cols-2*r-1, ...]    imDst[:, cols-r: cols, ...] = np.tile(imCum[:, cols-1:cols, ...], tile) - imCum[:, cols-2*r-1 : cols-r-1, ...]
    return imDst
def _gf_color(I, p, r, eps, s=None):    """ Color guided filter    I - guide image (rgb)    p - filtering input (single channel)    r - window radius    eps - regularization (roughly, variance of non-edge noise)    s - subsampling factor for fast guided filter    """    fullI = I    fullP = p    if s is not None:        I = sp.ndimage.zoom(fullI, [1/s, 1/s, 1], order=1)        p = sp.ndimage.zoom(fullP, [1/s, 1/s], order=1)        r = round(r / s)
    h, w = p.shape[:2]    N = box(np.ones((h, w)), r)
    mI_r = box(I[:,:,0], r) / N    mI_g = box(I[:,:,1], r) / N    mI_b = box(I[:,:,2], r) / N
    mP = box(p, r) / N
    # mean of I * p    mIp_r = box(I[:,:,0]*p, r) / N    mIp_g = box(I[:,:,1]*p, r) / N    mIp_b = box(I[:,:,2]*p, r) / N
    # per-patch covariance of (I, p)    covIp_r = mIp_r - mI_r * mP    covIp_g = mIp_g - mI_g * mP    covIp_b = mIp_b - mI_b * mP
    # symmetric covariance matrix of I in each patch:    #       rr rg rb    #       rg gg gb    #       rb gb bb    var_I_rr = box(I[:,:,0] * I[:,:,0], r) / N - mI_r * mI_r;    var_I_rg = box(I[:,:,0] * I[:,:,1], r) / N - mI_r * mI_g;    var_I_rb = box(I[:,:,0] * I[:,:,2], r) / N - mI_r * mI_b;
    var_I_gg = box(I[:,:,1] * I[:,:,1], r) / N - mI_g * mI_g;    var_I_gb = box(I[:,:,1] * I[:,:,2], r) / N - mI_g * mI_b;
    var_I_bb = box(I[:,:,2] * I[:,:,2], r) / N - mI_b * mI_b;
    a = np.zeros((h, w, 3))    for i in range(h):        for j in range(w):            sig = np.array([                [var_I_rr[i,j], var_I_rg[i,j], var_I_rb[i,j]],                [var_I_rg[i,j], var_I_gg[i,j], var_I_gb[i,j]],                [var_I_rb[i,j], var_I_gb[i,j], var_I_bb[i,j]]            ])            covIp = np.array([covIp_r[i,j], covIp_g[i,j], covIp_b[i,j]])            a[i,j,:] = np.linalg.solve(sig + eps * np.eye(3), covIp)
    b = mP - a[:,:,0] * mI_r - a[:,:,1] * mI_g - a[:,:,2] * mI_b
    meanA = box(a, r) / N[...,np.newaxis]    meanB = box(b, r) / N
    if s is not None:        meanA = sp.ndimage.zoom(meanA, [s, s, 1], order=1)        meanB = sp.ndimage.zoom(meanB, [s, s], order=1)
    q = np.sum(meanA * fullI, axis=2) + meanB
    return q

def _gf_gray(I, p, r, eps, s=None):    """ grayscale (fast) guided filter        I - guide image (1 channel)        p - filter input (1 channel)        r - window raidus        eps - regularization (roughly, allowable variance of non-edge noise)        s - subsampling factor for fast guided filter    """    if s is not None:        Isub = sp.ndimage.zoom(I, 1/s, order=1)        Psub = sp.ndimage.zoom(p, 1/s, order=1)        r = round(r / s)    else:        Isub = I        Psub = p

    (rows, cols) = Isub.shape
    N = box(np.ones([rows, cols]), r)
    meanI = box(Isub, r) / N    meanP = box(Psub, r) / N    corrI = box(Isub * Isub, r) / N    corrIp = box(Isub * Psub, r) / N    varI = corrI - meanI * meanI    covIp = corrIp - meanI * meanP

    a = covIp / (varI + eps)    b = meanP - a * meanI
    meanA = box(a, r) / N    meanB = box(b, r) / N
    if s is not None:        meanA = sp.ndimage.zoom(meanA, s, order=1)        meanB = sp.ndimage.zoom(meanB, s, order=1)
    q = meanA * I + meanB    return q

def _gf_colorgray(I, p, r, eps, s=None):    """ automatically choose color or gray guided filter based on I's shape """    if I.ndim == 2 or I.shape[2] == 1:        return _gf_gray(I, p, r, eps, s)    elif I.ndim == 3 and I.shape[2] == 3:        return _gf_color(I, p, r, eps, s)    else:        print("Invalid guide dimensions:", I.shape)

def guided_filter(I, p, r, eps, s=None):    """ run a guided filter per-channel on filtering input p        I - guide image (1 or 3 channel)        p - filter input (n channel)        r - window raidus        eps - regularization (roughly, allowable variance of non-edge noise)        s - subsampling factor for fast guided filter    """    if p.ndim == 2:        p3 = p[:,:,np.newaxis]
    out = np.zeros_like(p3)    for ch in range(p3.shape[2]):        out[:,:,ch] = _gf_colorgray(I, p3[:,:,ch], r, eps, s)    return np.squeeze(out) if p.ndim == 2 else out

def test_gf():    import imageio    cat = imageio.imread('cat.bmp').astype(np.float32) / 255    tulips = imageio.imread('tulips.bmp').astype(np.float32) / 255
    r = 8    eps = 0.05
    cat_smoothed = guided_filter(cat, cat, r, eps)    cat_detail = cat / cat_smoothed    print(cat_detail.shape)    cat_smoothed_s4 = guided_filter(cat, cat, r, eps, s=4)
    imageio.imwrite('cat_smoothed.png', cat_smoothed)    imageio.imwrite('cat_smoothed_s4.png', cat_smoothed_s4)    imageio.imwrite('cat_smoothed_detailed.png',cat_detail)
    tulips_smoothed4s = np.zeros_like(tulips)    tulips_detailed = np.zeros_like(tulips)    for i in range(3):        tulips_smoothed4s[:,:,i] = guided_filter(tulips, tulips[:,:,i], r, eps, s=4)
    tulips_detailed = tulips / tulips_smoothed4s    imageio.imwrite('tulips_detailed.png',tulips_detailed)    imageio.imwrite('tulips_smoothed4s.png', tulips_smoothed4s)
    tulips_smoothed = np.zeros_like(tulips)    for i in range(3):        tulips_smoothed[:,:,i] = guided_filter(tulips, tulips[:,:,i], r, eps)    imageio.imwrite('tulips_smoothed.png', tulips_smoothed)
if __name__ == '__main__':    test_gf()

一副圖像，經(jīng)過mask是圖像本身的導(dǎo)向濾波之后，得到一張細節(jié)圖和一張濾波圖。下面從左到右分別是原圖，細節(jié)圖和濾波圖。其實這是現(xiàn)在很多l(xiāng)ow-level領(lǐng)域的預(yù)處理步驟。拿到細節(jié)圖之后可以用卷積神經(jīng)網(wǎng)絡(luò)做下面的處理。

這里還推薦一個很好的輪子，C++實現(xiàn)的。

https://github.com/atilimcetin/guided-filter

結(jié)語 

以上就是常見的四種濾波算法的介紹。