目 錄
-
-
-
- 一、基本概念及原理
-
- 1. 掩模 Mask
- 2. 中值濾波 Median Filtering
- 3. 高斯濾波 Gaussian Filtering
- 4. 反銳化掩模 Unsharp Masking
- 5. Sobel邊緣檢測
- 6.Laplacian 邊緣檢測
- 7.Canny算法
- 二、具體實作及結果（MATLAB）
-
- 第 1 題
-
- 1.中值濾波
- 2. 高斯濾波
- 3. 分析比較
- 第 2 題
-
- 1. Unsharp Masking
- 2. Sobel Edge Detector
- 3. Laplacian Edge Detector
- 4. Canny Algorithm
- 參考資料
-
-
一、基本概念及原理
1. 掩模 Mask
掩模是數字圖像濾波中最基本的操作，就是用一個視窗在圖像上進行滑動，并對視窗内的像素點進行各種運算，這個視窗就是掩模（mask）。常用的掩模是d×d的正方形（為保證有中心點，d一般為奇數）。
2. 中值濾波 Median Filtering
中值濾波是一種非線性低通濾波。所謂中值濾波，就是對于每個像素點，找到其掩模覆寫範圍内像素點灰階值的中位數，作為這一點的灰階值。這樣處理可以使得圖像更加平滑，且效果較為良好，是實際中廣泛使用的一種平滑濾波方式。
3. 高斯濾波 Gaussian Filtering
高斯濾波是一種線性低通濾波，相當于在每個像素點的掩模中，按照二維高斯分布的方式增加了權重，将掩模中的點按照權重相加即得到中心點的值。
二維高斯機率密度函數如下：
f ( x , y ) = 1 2 π σ 2 e − x 2 + y 2 2 σ 2 f(x,y)=\frac{1}{2\pi\sigma^2}e^{-\frac{x^2+y^2}{2\sigma^2}} f(x,y)=2πσ21e−2σ2x2+y2
将掩模中心點坐标設為（0,0），其他點坐标分别在此基礎上加減即可得到，可以看出掩模中所有點坐标相加為0
由此得到掩模中每個位置的權重：
w ( i , j ) = 1 2 π σ 2 e − ( i − d − 1 ) 2 + ( j − d − 1 ) 2 2 σ 2 w(i,j)=\frac{1}{2\pi\sigma^2}e^{-\frac{(i-d-1)^2+(j-d-1)^2}{2\sigma^2}} w(i,j)=2πσ21e−2σ2(i−d−1)2+(j−d−1)2
其中 d d d是掩模半徑（掩模大小為 ( 2 d + 1 ) × ( 2 d + 1 ) (2d+1)×(2d+1) (2d+1)×(2d+1)），由此就可以得到高斯濾波模闆 W W W
根據機率密度的性質，這樣計算得到的權重之和為1，恰好能使得對圖像的灰階處理不超出限度且變換後的總體灰階等級沒有改變。
數字濾波中常使用二維卷積來實作對每個點進行模闆中權重去和的操作，此處即：
A ′ = A ∗ W A'=A*W A′=A∗W
A A A、 A ′ A' A′分别是濾波前和濾波後的圖像灰階矩陣。
4. 反銳化掩模 Unsharp Masking
反銳化掩模是先将圖像進行低通濾波（常用高斯濾波），然後将原圖像與之作差得到保留高頻成分的圖像，再将這個圖像與原圖像相加，得到增強邊緣的圖像（因為圖像邊緣往往灰階值變化率較快，即頻率較大）。反銳化掩模可以起到提高圖像高頻成分，增強圖像邊緣輪廓的作用。
其中，得到的保留高頻成分的圖像有時常乘以一個系數 k 再與原圖像相加，當 k = 1 k=1 k=1 時就是反銳化掩模，當 k ＞ 1 k＞1 k＞1時，稱為高提升濾波（high boost filtering）。
5. Sobel邊緣檢測
通過Sobel算子與原圖像卷積即可。通過這樣的操作可以得到圖像的邊緣，且這種方法得到的邊緣往往較為粗大。Sobel算子用于得到圖像的梯度資訊，則這樣處理後梯度較大的位置更亮，而梯度較大的位置往往是圖像的邊緣位置，故而通過這種方法可以實作圖像的邊緣檢測和提取。
Sobel算子可以得到圖像的梯度資訊，分為兩個部分：
x方向Sobel算子： [ 1 0 − 1 2 0 − 2 1 0 − 1 ] \begin{bmatrix}1&0&-1\\2&0&-2\\1&0&-1\end{bmatrix} ⎣⎡121000−1−2−1⎦⎤
y方向Sobel算子： [ 1 2 1 0 0 0 − 1 − 2 − 1 ] \begin{bmatrix}1&2&1\\0&0&0\\-1&-2&-1\end{bmatrix} ⎣⎡10−120−210−1⎦⎤
将它們分别與圖像卷積得到x，y方向上的梯度矩陣：
G x = [ 1 0 − 1 2 0 − 2 1 0 − 1 ] ∗ A G_x=\begin{bmatrix}1&0&-1\\2&0&-2\\1&0&-1\end{bmatrix}*A Gx=⎣⎡121000−1−2−1⎦⎤∗A ； G y = [ 1 2 1 0 0 0 − 1 − 2 − 1 ] ∗ A G_y=\begin{bmatrix}1&2&1\\0&0&0\\-1&-2&-1\end{bmatrix}*A Gy=⎣⎡10−120−210−1⎦⎤∗A
綜合兩個方向後有： G = G x 2 + G y 2 G=\sqrt{G_x^2+G_y^2} G=Gx2+Gy2
常将上式簡化為： ∣ G ∣ = ∣ G x ∣ + ∣ G y ∣ |G|=|G_x|+|G_y| ∣G∣=∣Gx∣+∣Gy∣
故而梯度方向矩陣也就容易求出了： Θ = a r c t a n ( G y G x ) \Theta=arctan\left(\frac{G_y}{G_x}\right) Θ=arctan(GxGy)
6.Laplacian 邊緣檢測
一般來說，圖像的邊緣輪廓不僅具有較高的梯度，也具有較高的二階導數。對一個二維函數求二階導的變換就是Laplacian變換，是以我們容易想到通過構造一個Laplace算子來實作圖像的邊緣檢測，即Laplacian邊緣檢測。
下面給出Laplace算子的推導過程：
對于一個二維函數，其Laplacian變換為： f ( x , y ) = ∂ 2 f ∂ x 2 + ∂ 2 f ∂ y 2 f(x,y)=\frac{\partial^2f}{\partial{x^2}}+\frac{\partial^2f}{\partial{y^2}} f(x,y)=∂x2∂2f+∂y2∂2f
其中
∂ 2 f ∂ x 2 = ∂ ∂ x [ f ( x + △ x 1 , y ) − f ( x , y ) ] △ x 1 2 = ∂ ∂ x [ f ( x + 1 , y ) − f ( x , y ) ] = [ f ( x + 1 + △ x 2 , y ) − f ( x + △ x 2 , y ) ] − [ f ( x + 1 , y ) − f ( x , y ) ] △ x 2 = f ( x , y ) − f ( x − 1 , y ) − f ( x + 1 , y ) + f ( x , y ) − 1 = f ( x + 1 , y ) + f ( x − 1 , y ) − 2 f ( x , y ) \frac{\partial^2f}{\partial{x^2}}=\frac{\partial}{\partial{x}}\frac{[f(x+\triangle{x_1},y)-f(x,y)]}{\triangle{x_1^2}} \\ \\=\frac{\partial}{\partial{x}}[f(x+1,y)-f(x,y)] \\ \\=\frac{[f(x+1+\triangle{x_2},y)-f(x+\triangle{x_2},y)]-[f(x+1,y)-f(x,y)]}{\triangle{x_2}} \\ \\=\frac{f(x,y)-f(x-1,y)-f(x+1,y)+f(x,y)}{-1} \\ \\=f(x+1,y)+f(x-1,y)-2f(x,y) ∂x2∂2f=∂x∂△x12[f(x+△x1,y)−f(x,y)] =∂x∂[f(x+1,y)−f(x,y)] =△x2[f(x+1+△x2,y)−f(x+△x2,y)]−[f(x+1,y)−f(x,y)] =−1f(x,y)−f(x−1,y)−f(x+1,y)+f(x,y) =f(x+1,y)+f(x−1,y)−2f(x,y)
因為數字圖像都是離散值，故 △ x 1 \triangle{x_1} △x1和 △ x 2 \triangle{x_2} △x2分别取 1 和 -1 （實際上就是二階差分）
同理， ∂ 2 f ∂ y 2 = f ( x , y + 1 ) + f ( x , y − 1 ) − 2 f ( x , y ) \frac{\partial^2f}{\partial{y^2}}=f(x,y+1)+f(x,y-1)-2f(x,y) ∂y2∂2f=f(x,y+1)+f(x,y−1)−2f(x,y)
是以，對于數字圖像這樣離散的矩陣來說，其Laplacian變換近似為：
∇ 2 f = [ f ( x + 1 , y ) + f ( x − 1 , y ) + f ( x , y + 1 ) + f ( x , y − 1 ) ] − 4 f ( x , y ) \nabla^2f=[f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)]-4f(x,y) ∇2f=[f(x+1,y)+f(x−1,y)+f(x,y+1)+f(x,y−1)]−4f(x,y)
為了友善做卷積，得到矩陣形式的Laplace算子：
[ 0 1 0 1 − 4 1 0 1 0 ] \begin{bmatrix}0&1&0\\1&-4&1\\0&1&0\end{bmatrix} ⎣⎡0101−41010⎦⎤
将其在8-鄰域進行擴充，得到擴充的Laplace算子：
[ 1 1 1 1 − 8 1 1 1 1 ] \begin{bmatrix}1&1&1\\1&-8&1\\1&1&1\end{bmatrix} ⎣⎡1111−81111⎦⎤
這即是通常使用的Laplace模闆
7.Canny算法
Canny算法是對Sobel算法的更進一步操作。首先用高斯濾波降低噪音對圖像的影響，再用Sobel算子計算出圖像的梯度，包括幅值和方向。對每個像素點，判斷其是否是其梯度方向上的極大值，如果不是則抑制（置0）。最後用雙門檻值進行檢測和連接配接，即大于高門檻值的确定保留，小于低門檻值的一定抑制，介于兩者之間的像素點，如果它的8-鄰域中存在确定保留點，則保留，否則抑制。
其中判斷梯度方向極大值時，需要使用插值法建構虛拟像素點，稱為亞像素點。使用Canny算法得到的圖像邊緣很細，且對噪音的濾去效果較好。
二、具體實作及結果（MATLAB）
第 1 題
空域低通濾波器：分别用高斯濾波器和中值濾波器去平滑測試圖像test1和2，模闆大小分别是3x3 ， 5x5 ，7x7； 分析各自優缺點
注：利用固定方差 sigma=1.5産生高斯濾波器.
1.中值濾波
實作步驟：
- 設定掩模大小，并使用
padarray()
- 對每個像素點，在掩模範圍内找灰階值的中位數賦給它
函數代碼：
function Img_out=MedianFilter(Img,masksize)
exsize=floor(masksize/2); %各方向擴充大小
Imgex=padarray(Img,[exsize,exsize],'replicate','both'); %擴充圖檔
[m,n]=size(Img);
Img_out=Img; %将Img_out準備為和Img相同的size
for i=1:m
for j=1:n
neighbor=Imgex(i:i+masksize-1,j:j+masksize-1); %截取鄰域
Img_out(i,j)=median(neighbor(:)); %中值濾波
end
end
end
得到結果：
2. 高斯濾波
實作步驟：
- 以掩模中心點為原點，得到掩模内各點的坐标
- 按照前述公式，計算出每個點的高斯權重，得到高斯模闆
- 将原圖與高斯模闆進行二維卷積得到處理後的圖像
函數代碼：
function Img_out=GaussianFilter(Img,masksize,sigma)
for i=1:masksize
for j=1:masksize
x=i-ceil(masksize/2);
y=j-ceil(masksize/2);
h(i,j)=exp(-(x^2+y^2)/(2*sigma^2))/(2*pi*sigma^2);
end
end
Img_out=uint8(conv2(Img,h,'same'));
end
得到結果：
3. 分析比較
中值濾波的優點是計算簡單迅速，能夠有效的去除圖像中的椒鹽噪聲；缺點則是對于一些高斯型噪聲不能很好的濾除。
而高斯濾波的優點則在于對高斯型噪聲有非常好的降噪效果，缺點則是計算相對複雜。同時，高斯濾波相當于對整幅圖像進行了模糊化處理，這一點有利有弊，視情況而定。
對于這幾幅圖來看，采用中值濾波的效果略優于高斯濾波的效果。這兩種濾波效果都與掩模大小密切相關，掩模越大，濾波效果越好，而且圖像更加明亮
第 2 題
利用高通濾波器濾波測試圖像test3,4：包括unsharp masking, Sobel edge detector, Laplace edge detection；Canny algorithm.分析各自優缺點
1. Unsharp Masking
實作步驟：
- 用高斯濾波得到模糊化圖像 Imgblur
- 原圖像Img-Imgblur 得到反銳化掩模 UnsharpMask
- UnsharpMask + Img 得到反銳化掩模處理後的圖像
- k*UnsharpMask + Img 得到高提升濾波後的圖像（這裡選 k = 3）
函數代碼：
function [Img_unsharp,Img_highboost,mask_unsharp]=Unsharp(Img)
Imgblur=GaussianFilter(Img,7,3); %用高斯濾波得到虛化圖像
mask_unsharp=Img-Imgblur; %得到反銳化掩模
Img_unsharp=Img+mask_unsharp; %得到反銳化掩模處理後的圖像
Img_highboost=Img+3*mask_unsharp; %得到highboost濾波後的圖像（k=3）
end
得到結果：
算法分析：
如原理部分所述，反銳化掩模、高提升濾波算法可以使得圖像的邊緣輪廓和高頻成分增強，其效果與 k 取值密切相關，k 越大對細節的增強就越明顯。
2. Sobel Edge Detector
實作步驟：
- 構造x，y方向的Sobel算子Sx和Sy
- 分别用原圖與Sx，Sy進行二維卷積，得到兩個方向的梯度矩陣
- 将兩個方向梯度分别求絕對值并相加得到圖像的梯度矩陣
- 将梯度矩陣轉成uint8格式，即為處理後圖像
函數代碼：
function [Img_out,gx,gy]=Sobel(Img)
sobely=[1,2,1;0,0,0;-1,-2,-1]; %構造y方向Sobel算子
sobelx=[1,0,-1;2,0,-2;1,0,-1]; %構造x方向Sobel算子
gx=conv2(Img,sobelx,'same');
gy=conv2(Img,sobely,'same');
Img_out=uint8(abs(gx)+abs(gy));
end
得到結果：
算法分析：
采用Sobel算法提取出的邊緣範圍較大，是以顯示的邊緣較為粗大，看起來更加清晰明了，且對噪聲有不錯的抑制效果，但是精确度不高。
3. Laplacian Edge Detector
實作步驟：
- 構造出擴充Laplace算子
- 将原圖與Laplace算子進行二維卷積，結果轉為uint8格式即為處理後圖像
函數代碼：
function Img_out=Laplace(Img)
laplace=[1,1,1;1,-8,1;1,1,1]; %構造拉普拉斯算子
Img_out=uint8(conv2(Img,laplace,'same'));
end
得到結果：
算法分析：
采用Laplacian邊緣檢測算法，可以看出當圖像比較"幹淨"時有較好的邊緣檢測作用，且精度明顯優于Sobel算法。但是Laplacian對噪聲過于敏感，很容易受到噪音幹擾，可以看到即使test3的檢測結果仍然有很多椒鹽噪聲。
4. Canny Algorithm
實作步驟：
- 用高斯濾波将圖像模糊化處理以降低噪音影響
- 用Sobel算子得到x，y方向的梯度幅值和方向
- 對梯度方向非極大值進行抑制
- 将梯度方向從0°開始，每45°化為一個區域，共四個區域（正負視為同一方向）
- 對像素點所在區域找到區域兩側的正像素點，根據梯度方向進行權重平均，就可以插值得到兩個亞像素點g1，g2，分别處于梯度方向的兩側
- 判斷像素點是否大于g1，g2，即判斷在梯度方向上是否為極大值，若不是則置0
- 設定高低門檻值進行雙門檻值檢測
- 大于高門檻值T2的值一定保留，小于低門檻值T1的一定置0
- 介于中間的點，如果在其8-鄰域内有确定保留的點，則它保留，否則置0
- 由于Canny得到的邊緣往往較細，可以對整幅圖進行強化處理
函數代碼：
function Img_out=Canny(Img)
Imgblur=GaussianFilter(Img,3,2); %用高斯濾波對圖像進行模糊處理以降低噪聲影響
[~,gx,gy]=Sobel(Imgblur); %用Sobel得到x，y方向的梯度
g=uint8(sqrt(double(gx.*gx)+double(gy.*gy))); %得到圖像的梯度幅值矩陣
theta=atan2(double(gy),double(gx)); %得到圖像的梯度方向矩陣
gcom=g; %準備一個存儲空間gcom
g=double(g); %将g化為double型，為下面運算準備
%下面将按照方向分别判斷：
for i=2:size(g,1)-1
for j=2:size(g,2)-2
if ((theta(i,j)>=0)&(theta(i,j)<pi/4))|((theta(i,j)>=-pi)&(theta(i,j)<-3*pi/4))
g1=g(i-1,j+1)*(gy(i,j)/gx(i,j))+g(i,j+1)*(1-gy(i,j)/gx(i,j)); %用插值得到梯度方向上的一個亞像素點，下同
g2=g(i+1,j-1)*(gy(i,j)/gx(i,j))+g(i,j-1)*(1-gy(i,j)/gx(i,j)); %用插值得到梯度方向上的另一個亞像素點，下同
if (gcom(i,j)<=g1)|(gcom(i,j)<=g2)
gcom(i,j)=0; %若在梯度方向上不是極大值則置0，下同
end
end
if ((theta(i,j)>=pi/4)&(theta(i,j)<pi/2))|((theta(i,j)>=-3*pi/4)&(theta(i,j)<-pi/2))
g1=g(i-1,j+1)*(gx(i,j)/gy(i,j))+g(i-1,j)*(1-gx(i,j)/gy(i,j));
g2=g(i+1,j-1)*(gx(i,j)/gy(i,j))+g(i+1,j)*(1-gx(i,j)/gy(i,j));
if (gcom(i,j)<=g1)|(gcom(i,j)<=g2)
gcom(i,j)=0;
end
end
if ((theta(i,j)>=pi/2)&(theta(i,j)<3*pi/4))|((theta(i,j)>=-pi/2)&(theta(i,j)<-pi/4))
g1=g(i-1,j-1)*(-gx(i,j)/gy(i,j))+g(i-1,j)*(1+gx(i,j)/gy(i,j));
g2=g(i+1,j+1)*(-gx(i,j)/gy(i,j))+g(i+1,j)*(1+gx(i,j)/gy(i,j));
if (gcom(i,j)<=g1)|(gcom(i,j)<=g2)
gcom(i,j)=0;
end
end
if ((theta(i,j)>=3*pi/4)&(theta(i,j)<=pi))|((theta(i,j)>=-pi/4)&(theta(i,j)<0))
g1=g(i-1,j-1)*(-gy(i,j)/gx(i,j))+g(i,j-1)*(1+gy(i,j)/gx(i,j));
g2=g(i+1,j+1)*(-gy(i,j)/gx(i,j))+g(i,j+1)*(1+gy(i,j)/gx(i,j));
if (gcom(i,j)<=g1)|(gcom(i,j)<=g2)
gcom(i,j)=0;
end
end
end
end
g=uint8(gcom); %将更新後的矩陣存回g中
high=30;low=5; %設定高低門檻值
gunderlow=uint8(g>=low); %低于低門檻值的點
g1=g.*gunderlow; %将低于低門檻值的點置0
gbetween=uint8(g1<=high); %介于兩門檻值之間的點（待确定點）
gsure=uint8(g1>high); %高于高門檻值的點，即确定保留的點
[m,n]=size(gsure);z1=zeros(1,n);z2=zeros(m,1); %z1，z2分别是水準豎直的0向量
%以下操作将高于高門檻值的點擴充到8-鄰域：
gsureex=gsure+[z2,gsure(:,1:n-1)]+[gsure(:,2:n),z2]+[z1;gsure(1:m-1,:)]+[gsure(2:m,:);z1]...
+[z2,[z1(1:n-1);gsure(1:m-1,1:n-1)]]+[z2,[gsure(2:m,1:n-1);z1(1:n-1)]]...
+[[z1(1:n-1);gsure(1:m-1,2:n)],z2]+[[gsure(2:m,2:n);z1(1:n-1)],z2];
gsureex=uint8(gsureex>0);
gbetween1=gbetween.*gsureex; %得到8-鄰域中有高于高門檻值點的待确定點
Img_out=g1.*(gbetween1+gsure); %兩者合并即為最終保留的點
Img_out=255*Img_out; %将最終保留的點強化
end
得到結果：
算法分析：
Canny算法幾乎不受噪音幹擾，且得到的邊緣精度很高，但是算法較為複雜費時，并且處理效果與高斯濾波的兩個參數（掩模大小、标準差sigma），以及自身的兩個參數（高低門檻值）的選取密切相關，這幾個參數的選取也根據圖像有所不同。在本題中，為将兩幅圖一起處理輸出，綜合比較後選取高斯掩模大小為3，sigma為2，高門檻值為30，低門檻值為5，得到的效果較為良好。
參考資料
- 岡薩雷斯.數字圖像處理（第三版）[M].北京:電子工業出版社.2017.1
數字圖像處理作業連結：
數字圖像處理第一次作業——Bmp格式與基本變換
數字圖像處理第二次作業——圖像仿射變換
數字圖像處理第三次作業——基于直方圖的圖像空域操作
數字圖像處理第四次作業——空域濾波與邊緣檢測
數字圖像處理第五次作業——頻域濾波器
數字圖像處理第六次作業——圖像噪聲和恢複
數字圖像處理第七次作業——圖像直線檢測