詳解數字圖像的濾波和邊緣檢測

一、圖像濾波

圖像濾波的主要目的就是在盡量保留圖像細節特征的條件下對目標圖像的噪聲進行抑制。圖像濾波其主要分為線性濾波器和非線性濾波器。

1.1 線性濾波器

1.1.1 均值濾波

原理：在圖像上，對待處理的像素給定一個模板，該模板包括了其周圍的鄰近像素。將模板中的全體像素的均值來替代原來的像素值的方法。以 \(3\times3\) 的均值濾波器的模板為：

\[w= \begin{bmatrix} 1&1&1\\ 1&1&1\\ 1&1&1 \end{bmatrix} \]

模板所覆蓋的圖像的像素為：

\[Image = \begin{bmatrix} (i-1,j-1)&(i-1,j)&(i-1,j+1)\\ (i,j-1)&(i,j)&(i,j+1)\\ (i+1,j-1)&(i+1,j-1)&(i+1,j+1) \end{bmatrix} \]

因此，模板中心點所對應的像素值可以表示為：

\[g(i,j) = \dfrac{1}{N} \sum_{i=1}^{N}Image \odot w = \dfrac{1}{N}\sum_{i=1}^{N}\sum_{m=-1}^{1}\sum_{n=-1}^{1}f(i+m,j+n) \]

在進行濾波之前，一般會先對原圖像進行填充，填充的像素大小為：\(r = floor(模板寬度/2)\)，不然的話始終無法對圖像邊緣的像素進行濾波。

均值濾波的特點：均值濾波能夠降低圖像中的尖銳變化，去除圖像中的噪聲，且濾波窗口越大，去噪效果越好同時圖像也會變得越模糊

1.1.1 高斯濾波

原理：其實高斯濾波和均值濾波的原理差不多，都是用一個模板（或稱卷積、掩模）掃描圖像中的每一個像素，用模板確定的鄰域內像素的加權平均灰度值去替代模板中心像素點的值。區別就在於，高斯濾波中模板的數值是通過高斯函數生成的。還是以 \(3\times3\) 的模板為例，以模板的左上角為坐標原點，模板中心的坐標為 \((1,1)\)，則生成的高斯模板為：

\[w(i,j) = \dfrac{1}{2\pi \sigma^{2}}exp{-\dfrac{(i-1)^2+(j-1)^2}{2\sigma^{2}}} \]

因此，模板中心點所對應的像素值可以表示為：

\[g(i,j) = Image \odot \begin{bmatrix} w_{(i-1,j-1)}&w_{(i-1,j)}&w_{(i-1,j+1)}\\ w_{(i,j-1)}&w_{(i,j)}&w_{(i,j+1)}\\ w_{(i+1,j-1)}&w_{(i+1,j)}&w_{(i+1,j+1)} \end{bmatrix} = \sum_{m=-1}^{1}\sum_{n=-1}^{1} w(i+m,j+n)f(i+m,j+n) \]

高斯濾波的特點：根據公式可以得到，在計算高斯濾波模板時，高斯濾波模板的值由中間向四周遞減，且標准差越小，二維高斯圖像越窄小，平滑效果不明顯；標准差越大，而為高斯圖像越矮寬，濾波效果比較明顯。下圖展示了不同 \(\sigma\) 的\(7\times7\) 高斯模板的差別：

1.2 非線性濾波器

1.2.1 中值濾波器

原理：中值濾波器就是取模板覆蓋區域的排序之后的中間值作為該模板區域內中心的像素值。還是以 \(3\times 3\) 中值濾波器模板為例，此時覆蓋的像素為：

\[Image = \begin{bmatrix} (i-1,j-1)&(i-1,j)&(i-1,j+1)\\ (i,j-1)&(i,j)&(i,j+1)\\ (i+1,j-1)&(i+1,j-1)&(i+1,j+1) \end{bmatrix} \]

那么，中值濾波后的像素值為：

\[g(i,j) = median(Image) \]

中值濾波器的特點：中值濾波器很適合處理椒鹽噪聲。

1.2.2 雙邊濾波器

設計原因：在均值濾波器和高斯濾波器時，濾波器的模板值僅僅考慮到了像素的空間位置，潛在的假設為：離模板中心點越近，那么他們的像素值也越接近。但是在紋理比較復雜的區域，這種假設很難成立，因此采用上述濾波器進行濾波時，圖像的邊緣信息保持不夠理想。那么能不能設計一款濾波器，能夠識別紋理復雜的區域，從而實現更好的邊緣保持i效果？

原理：雙邊濾波器分為兩個部分，首先采用高斯濾波器生成只跟像素位置相關的模板，然后根據像素值，生成另外一個跟像素值相關的模板。兩個模板逐像素點擊就形成了雙邊濾波模板。具體描述如下：

\[w1(i,j,k) = exp(-\dfrac{(i-k)^2+(j-k)^2}{2\sigma_d^{2}}) \\ w2(i,j,k) = exp(-\dfrac{||f(i,j)-f(k,k)||^{2}}{2\sigma_d^{2}}) \]

因此，最終的雙邊濾波模板為：

\[w = w1 * w2 = exp(-\dfrac{(i-k)^2+(j-k)^2}{2\sigma_d^{2}}-\dfrac{||f(i,j)-f(k,k)||^{2}}{2\sigma_d^{2}}) \]

在雙邊濾波器中多考慮了值域模板，當兩個位置像素越接近，那么 \(w2\) 對應位置的權重就越大，否則就會越小。因此通過 \(w2\) 我們在濾波器中保護了圖像的邊緣信息。其模板的圖像如下：

雙邊濾波器的特點：雙邊濾波器是一種可以保邊緣信息的去噪濾波器。

1.3 濾波器的 Matlab 代碼實現

clear all
clc
image = imread('../測試圖像/2.bmp');
[H,W,C] = size(image);  %記錄讀入圖像的長寬

H_win = 7;  %定義模板窗口的寬度
image_input = double(rgb2ycbcr(image)); %轉換成Ycbcr顏色空間
image_input = image_input(:,:,1); %取Y通道來進行處理

paddle = floor(H_win/2);    %計算需要填充的像素寬度
image_data = zeros(H+2*paddle, W+2*paddle); %0填充

Start = paddle+1;   %計算圖像的開始和結束位置
End = H+paddle;
image_data(Start:End, Start:End) = image_input(:,:);%將圖像放在填充圖像的中間

%生成高斯濾波模板
sigma_g = 3;
sigma_r = 6;
[x,y] = meshgrid(-paddle:paddle,-paddle:paddle);
k_gaussian1=exp(-(x.^2+y.^2)/(2*sigma_g^2));     %以距離作為自變量高斯濾波器

k_gaussian=k_gaussian1/sum(k_gaussian1(:));%歸一化
Map_X = 1:H_win; Map_Y = 1:H_win; %作3維圖的橫坐標
surf(Map_X,Map_Y,k_gaussian)    %繪制高斯濾波核的圖


for i=Start:1:End
    for j=Start:1:End
        image_window = image_data(i-paddle:i+paddle, j-paddle:j+paddle);
        
        Mean_image(i-paddle,j-paddle) = mean(mean(image_window));   %均值濾波
        Median_image(i-paddle,j-paddle) = median(image_window(:));  %中值濾波
        
        Gaussian_image(i-paddle, j-paddle) = sum(sum(image_window.*k_gaussian));  %高斯濾波
        
        double_w = exp(-(image_window-image_data(i,j)).^2/(2*sigma_r^2)); %以周圍和當前像素灰度差值作為自變量的高斯濾波器
        double_w = k_gaussian1 .* double_w;
        double_w = double_w/sum(double_w(:));
        
        Double_image(i-paddle, j-paddle) = sum(sum(image_data(i-paddle:i+paddle, j-paddle:j+paddle).*double_w));  %雙邊濾波
    end
end
figure, imshow(uint8(image_input)), title('原圖')
subplot(2,2,1), imshow(uint8(Mean_image)), title('5*5均值')
subplot(2,2,2), imshow(uint8(Median_image)), title('中值濾波')
subplot(2,2,3), imshow(uint8(Gaussian_image)), title('5*5,1.5 高斯')
subplot(2,2,4), imshow(uint8(Double_image)), title('5*5,1.5 雙邊')

代碼輸出結果如下：

二、圖像邊緣檢測

邊緣檢測是為了將其周圍像素灰度有階躍變化的像素檢測出來，這些像素組成的集合就是該圖像的邊緣。比較常用的邊緣檢測方法就是考察每個像素在某個領域內灰度的變化，然后利用邊緣臨近一階或二階方向導數變化規律檢測邊緣，即邊緣檢測局部算法

2.1 一階邊緣檢測算子

2.1.1 Sobel 算子

原理：Sobel 算子包括兩組 \(3\times3\) 的矩陣，左邊的表示垂直，右邊的表示水平。將它與圖像作平面卷積，即可分別得出垂直及水平的亮度差分近似值。

\[G_x = \begin{bmatrix} -1&0&1\\-2&0&2\\-1&0&1\\\end{bmatrix}, G_y = \begin{bmatrix} 1&2&1\\0&0&0\\-1&-2&-1\\\end{bmatrix} \]

假設 \(O_x\) 和 \(O_{y}\) 分別用來表示被模板 \(G_x\) 和 \(G_{y}\) 卷積后的結果，那么最終的輸出圖像邊緣可以表示為：

\[O = \sqrt{O_x^{2}+O_y^{2}} \]

Sobel算子的特點：Soble算子在水平和垂直兩個方向上求導，得到的是圖像在X方法與Y方向梯度圖像。但是比較敏感，容易受影響，要通過高斯模糊（平滑）來降噪。

2.1.2 Canny 算子

原理：Canny算子是由計算機科學家John F. Canny於1986年提出的一種邊緣檢測算子，是目前理論上相對最完善的一種邊緣檢測算法。主要包括以下幾個步驟：

1. 高斯濾波對圖像進行預處理，提前弱化噪聲，增加邊緣提取准確度
2. Sobel 算子進行像素梯度計算，得到梯度強度矩陣
3. 非極大值抑制：非極大值像素梯度抑制起到將邊緣“瘦身”的作用。其基本方法是將當前像素梯度強度與沿正負梯度方向上的相鄰像素的梯度強度進行比較，若其最大（即為極值），則保留該像素為邊緣點，若不是最大，則對其進行抑制，不將其作為邊緣點。
4. 滯后閾值處理：定義一個高閾值和一個低閾值。梯度強度低於低閾值的像素點被抑制，不作為邊緣點；高於高閾值的像素點被定義為強邊緣，保留為邊緣點；處於高低閾值之間的定義為弱邊緣，留待進一步處理。
5. 孤立弱邊緣抑制：通常而言，由真實邊緣引起的弱邊緣像素點將連接到強邊緣像素點，而噪聲響應則未連接。通過查看弱邊緣像素及其8個鄰域像素，可根據其與強邊緣的連接情況來進行判斷。一般，可定義只要其中鄰域像素其中一個為強邊緣像素點，則該弱邊緣就可以保留為強邊緣，即真實邊緣點。

Canny檢測的特點：Canny 算子對於濾波參數和高、低閾值的選取還是較為敏感，使得實際應用過程中需要反復調試。

2.2 二階邊緣檢測算子

2.2.1 拉普拉斯算子

原理：拉普拉斯算子是二階微分算子，實際上就是圖像梯度的散度：

\[\nabla^{2}f = \dfrac{\partial^2f}{\partial x^2}+\dfrac{\partial^2f}{\partial y^2} \]

然而：

\[\dfrac{\partial^2f}{\partial x^2}=f(x,y+1)-2f(x,y)+f(x,y-1) \\ \dfrac{\partial^2f}{\partial y^2}=f(x+1,y)-2f(x,y)+f(x-1,y) \]

所以，卷積的模板可以寫為：

\[\nabla^{2}f = \begin{bmatrix}0&1&0\\1&-4&1\\0&1&0\end{bmatrix}, \nabla^{2}f = \begin{bmatrix}1&1&1\\1&-8&1\\1&1&1\end{bmatrix} \]

拉普拉斯算子的特點：Laplacian算子法對噪聲比較敏感，所以很少用該算子檢測邊緣，而是用來判斷邊緣像素視為與圖像的明區還是暗區。Laplacian算子是各向同性的，能對任何走向的界線和線條進行銳化，無方向性。

2.3 邊緣檢測的 Matlab 代碼實現

clear all
clc
image = imread('../測試圖像/2.bmp');
[H,W,C] = size(image);  %記錄讀入圖像的長寬

H_win = 21;  %定義模板窗口的寬度
image_input = double(rgb2ycbcr(image)); %轉換成Ycbcr顏色空間
image_input = image_input(:,:,1); %取Y通道來進行處理

%定義sobel算子
G_x = [-1,0,1;-2,0,2;-1,0,1];
G_y = [1,2,1;0,0,0;-1,-2,-1];
%定義拉普拉斯算子
L = [-1,-1,-1;-1,8,-1;-1,-1,-1];
image_data = mat2gray(image_input); %圖像歸一化

for i=2:1:H-1
    for j=2:1:W-1
        image_win = image_data(i-1:i+1, j-1:j+1);
        Gx = sum(sum(image_win .* G_x));
        Gy = sum(sum(image_win .* G_y));
        Sobel_image(i, j) = sqrt(Gx^2 + Gy^2);
        Lapla_image(i, j) = image_data(i,j+1)+image_data(i,j-1)+image_data(i+1,j)+image_data(i-1,j)-4*(image_data(i,j));
    end
end
subplot(131), imshow(image_data), title('原圖')
subplot(132), imshow(Sobel_image), title('Canny邊緣檢測')
subplot(133), imshow(Lapla_image), title('拉普拉斯檢測')

代碼輸出結果如下：

2.4 邊緣檢測算法總結

1. Sobel算子檢測方法對灰度漸變和噪聲較多的圖像處理效果較好，sobel算子對邊緣定位不是很准確，圖像的邊緣不止一個像素；當對精度要求不是很高時，是一種較為常用的邊緣檢測方法。
2. Canny方法不容易受噪聲干擾，能夠檢測到真正的弱邊緣。優點在於，使用兩種不同的閾值分別檢測強邊緣和弱邊緣，並且當弱邊緣和強邊緣相連時，才將弱邊緣包含在輸出圖像中。
3. Laplacian算子法對噪聲比較敏感，所以很少用該算子檢測邊緣，而是用來判斷邊緣像素視為與圖像的明區還是暗區。拉普拉斯高斯算子是一種二階導數算子，將在邊緣處產生一個陡峭的零交叉, Laplacian算子是各向同性的，能對任何走向的界線和線條進行銳化，無方向性。這是拉普拉斯算子區別於其他算法的最大優點。

三、參考資料

https://zhuanlan.zhihu.com/p/355263110

https://www.cnblogs.com/wangguchangqing/p/6416401.html

https://www.cnblogs.com/yibeimingyue/p/10878514.html

https://zhuanlan.zhihu.com/p/99959996

https://www.cnblogs.com/yibeimingyue/p/10878514.html

https://www.jianshu.com/p/2a06c68f6c14