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

ZhiboZhao發表於2021-08-18

一、影像濾波

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

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年提出的一種邊緣檢測運算元，是目前理論上相對最完善的一種邊緣檢測演算法。主要包括以下幾個步驟：

高斯濾波對影像進行預處理，提前弱化噪聲，增加邊緣提取準確度
Sobel 運算元進行畫素梯度計算，得到梯度強度矩陣
非極大值抑制：非極大值畫素梯度抑制起到將邊緣“瘦身”的作用。其基本方法是將當前畫素梯度強度與沿正負梯度方向上的相鄰畫素的梯度強度進行比較，若其最大（即為極值），則保留該畫素為邊緣點，若不是最大，則對其進行抑制，不將其作為邊緣點。
滯後閾值處理：定義一個高閾值和一個低閾值。梯度強度低於低閾值的畫素點被抑制，不作為邊緣點；高於高閾值的畫素點被定義為強邊緣，保留為邊緣點；處於高低閾值之間的定義為弱邊緣，留待進一步處理。
孤立弱邊緣抑制：通常而言，由真實邊緣引起的弱邊緣畫素點將連線到強邊緣畫素點，而噪聲響應則未連線。通過檢視弱邊緣畫素及其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 邊緣檢測演算法總結

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