圖像邊緣檢測各種算子MATLAB實現以及實際應用

1、圖像處理中的數學方法實驗報告學生姓名：趙芳舟教師姓名：曾理學院：數學與統(tǒng)計學院專業(yè)：信息與計算科學學號：聯(lián)系方式：梯度和拉普拉斯算子在圖像邊緣檢測中的應用1、 數學方法邊緣檢測最通用的方法是檢測灰度值的不連續(xù)性，這種不連續(xù)性用一階和二階導數來檢測。1. （1）一階導數：一階導數即為梯度，對于平面上的圖像來說，我們只需用到二維函數的梯度，即： f=gxgy=fxfy,該向量的幅值：f=magf=gx2+gy21/2=f/x2+f/y21/2,為簡化計算，省略上式平方根，得到近似值fgx2+gy2；或通過取絕對值來近似，得到：fgx+gy。（2）二階導數：二階導數通常用拉普拉斯算子來計算，由二階微

2、分構成：2fx,y=2f(x,y)x2+2f(x,y)y22. 邊緣檢測的基本思想：（1） 尋找灰度的一階導數的幅度大于某個指定閾值的位置；（2） 尋找灰度的二階導數有零交叉的位置。3. 幾種方法簡介（1） Sobel邊緣檢測器：以差分來代替一階導數。Sobel邊緣檢測器使用一個33鄰域的行和列之間的離散差來計算梯度，其中，每行或每列的中心像素用2來加權，以提供平滑效果。f=gx2+gy21/2=z7+2z8+z9-(z1+2z2+z3)2+z3+2z6+z9-z1+2z4+z721/2-1-21000121 -101-202-101 （2） Prewitt邊緣檢測器:使用下圖所示模板來數字化

3、地近似一階導數。與Sobel檢測器相比，計算上簡單一些，但產生的結果中噪聲可能會稍微大一些。gx=z7+z8+z9-z1+z2+z3gy=z3+z6+z9-z1-z4-z7-1-1-1000111-101-101-101（3） Roberts邊緣檢測器：使用下圖所示模板來數字化地將一階導數近似為相鄰像素之間的差，它與前述檢測器相比功能有限（非對稱，且不能檢測多種45倍數的邊緣）。gx=z9-z5gy=z8-z6-10010-110（4） Laplace邊緣檢測器：二維函數f(x,y)的拉普拉斯是一個二階的微分定義：2fx,y=2f(x,y)x2+2f(x,y)y2模板算子可分為四鄰域和八鄰域，

4、如下：0101-41010（四鄰域）1111-81111（八鄰域）（5） LoG邊緣檢測器由于噪聲點（灰度與周圍點相差很大的像素點）對邊緣檢測有一定的影響，所以效果更好的是LoG算子，即Laplacian-Guass算子。引入高斯函數來平滑噪聲：該函數的Laplace算子：它把Guass平滑濾波器和Laplace銳化濾波器結合起來，先平滑掉噪聲，再進行邊緣檢測，所以效果比單用Laplace算子要更為平滑，效果更好。（6） Canny邊緣檢測器主要分為以下幾個步驟：使用具有指定標準差的一個高斯濾波器來平滑圖像，以減少噪聲；在每個點處計算局部梯度和邊緣方向；對步驟中確定的邊緣點產生梯度中的脊線頂部

5、進行追蹤，并將實際山不在脊線頂部的像素設置為零，從而在輸出中給出一條細線（非最大值抑制），然后使用滯后閾值處理法對這些脊線像素進行閾值處理。最后進行邊緣連接。2、 實驗結果原始圖像：1.Roberts算子2.Prewitt算子3.Sobel算子4.Laplace算子（4鄰域）（8鄰域）5. LoG四種結果的比較：（Laplace采用8鄰域模板）圖像來源：/show/.html（網行天下首頁設計圖庫文化藝術插畫集：牧羊犬灰度）圖像規(guī)格：800677 格式：SVG 模式：RGB硬件條件：軟件條件：運行環(huán)境：Matlab r2015b三、討論1. 各方法

6、優(yōu)缺點及比較：（1） Roberts算子：Roberts算子是22算子，對具有陡峭的低噪聲圖像響應最好，并且檢測垂直邊緣的效果好于斜向邊緣，定位精度高。然而，它對噪聲敏感，無法抑制噪聲的影響。因此，它適用于邊緣明顯且噪聲較少的圖像分割。（2） Prewitt算子：Prewitt算子將兩個點的各自一定領域內的灰度值求和，并根據兩個灰度值和的差來計算x,y的偏導數。它是平均濾波，對噪聲有抑制作用，對于灰度漸變的低噪聲圖像有較好的檢測效果，但是像素平均相當于對圖像的低通濾波，所以它對邊緣的定位不如Roberts算子。對于混合多復雜噪聲的圖像，效果不太理想。（3） Sobel算子：Sobel算子是濾波

7、算子的形式，用于提取邊緣，可以利用快速卷積函數， 簡單有效，因此應用廣泛。美中不足的是，Sobel算子并沒有將圖像的主體與背景嚴格地區(qū)分開來，即Sobel算子沒有嚴格地模擬人的視覺生理特征，所以提取的圖像輪廓有時并不能令人滿意。（4） Laplace算子：Laplace算子是一種各向同性算子，在只關心邊緣的位置而不考慮其周圍的象素灰度差值時比較合適。Laplace算子對孤立象素的響應要比對邊緣或線的響應要更強烈，因此只適用于無噪聲圖象。存在噪聲情況下，使用Laplacian算子檢測邊緣之前需要先進行低通濾波。（5） LoG算子：由于二階導數是線性運算，利用LoG卷積一幅圖像與首先使用高斯型平滑

8、函數卷積改圖像，然后計算所得結果的拉普拉斯是一樣的。所以在LoG公式中使用高斯函數的目的就是對圖像進行平滑處理，使用Laplacian算子的目的是提供一幅用零交叉確定邊緣位置的圖像；圖像的平滑處理減少了噪聲的影響并且它的主要作用還是抵消由Laplacian算子的二階導數引起的逐漸增加的噪聲影響。（6） Canny算子：是一個具有濾波，增強和檢測的多階段的優(yōu)化算子，在進行處理前，Canny算子先利用高斯平滑濾波器來平滑圖像以消除噪聲，Canny分割算法采用一階偏導的有限差分來計算梯度的幅值和方向。在處理過程中，Canny算法還將經過一個非極大值抑制的過程。最后Canny算法將采用兩個閾值來連接邊

9、緣。高定位精度、低誤判率、抑制虛假邊緣，適用于高噪聲圖像。四、實際應用1、概述：圖像邊緣檢測廣泛應用于車牌識別，人臉識別，地震帶檢測，生物醫(yī)療以及產品外觀檢測等方面。為了體現其實際應用，在此針對車牌識別來進行實驗，實驗流程如下：確定區(qū)域文字分割模板配對均值濾波邊緣檢測二值化灰度處理3、 實驗結果原圖像：灰度處理：將彩色圖像轉化為灰度圖像，縮小圖片占用的空間，減少處理時間邊緣檢測：車牌區(qū)域的邊緣信息一般比其他部分豐富且集中，只要 背景不太復雜，沒有過多干擾，邊緣特征可將車牌區(qū)域與其他部分區(qū)別開來，從而消除圖像中的無關信息。腐蝕和膨脹：腐蝕可以分割獨立的圖像元素，膨脹用于連接相鄰的元素，更加利于圖

10、像的分割最后將車牌號分割成七個單個字符，建立模板庫，與其進行配對，最終得出結果。五、程序介紹（1） Roberts邊緣檢測 clear;sourcePic=imread(C:Users34899Desktoptimg.jpg); %讀取原圖像grayPic=mat2gray(sourcePic); %實現圖像矩陣的歸一化操作m,n=size(grayPic);newGrayPic=grayPic;%為保留圖像的邊緣一個像素robertsNum=0; %經roberts算子計算得到的每個像素的值robertThreshold=0.2; %設定閾值for j=1:m-1 %進行邊界提取 for k

11、=1:n-1 robertsNum = abs(grayPic(j,k)-grayPic(j+1,k+1) + abs(grayPic(j+1,k)-grayPic(j,k+1); if(robertsNum robertThreshold) newGrayPic(j,k)=255; else newGrayPic(j,k)=0; end endendsubplot(1,2,1);imshow(sourcePic);title(原圖); subplot(1,2,2);imshow(newGrayPic);title(Robert算子處理后圖像); （2） Prewitt邊緣檢測 clear;s

12、ourcePic=imread(C:Users34899Desktoptimg.jpg); %讀取原圖像grayPic=mat2gray(sourcePic);%實現圖像矩陣的歸一化操作m,n=size(grayPic);newGrayPic=grayPic;%為保留圖像的邊緣一個像素PrewittNum=0;%經Prewitt算子計算得到的每個像素的值PrewittThreshold=0.5;%設定閾值for j=2:m-1 %進行邊界提取for k=2:n-1 PrewittNum=abs(grayPic(j-1,k+1)-grayPic(j+1,k+1)+grayPic(j-1,k)-g

13、rayPic(j+1,k)+grayPic(j-1,k-1)-grayPic(j+1,k-1)+abs(grayPic(j-1,k+1)+grayPic(j,k+1)+grayPic(j+1,k+1)-grayPic(j-1,k-1)-grayPic(j,k-1)-grayPic(j+1,k-1); if(PrewittNum PrewittThreshold) newGrayPic(j,k)=255; else newGrayPic(j,k)=0; end endendsubplot(1,2,1);imshow(sourcePic);title(原圖); subplot(1,2,2);ims

14、how(newGrayPic);title(Prewitt算子處理后圖像); （3） Sobel邊緣檢測 clear all; close all; imag = imread(C:Users34899Desktoptimg.jpg); %讀取關鍵幀 imag = rgb2gray(imag); %轉化為灰度圖 subplot(1,2,1);imshow(imag);title(原圖); high,width = size(imag); % 獲得圖像的高度和寬度 F2 = double(imag); U = double(imag); uSobel = imag; for i = 2:high

15、 - 1 %sobel邊緣檢測 for j = 2:width - 1 Gx = (U(i+1,j-1) + 2*U(i+1,j) + F2(i+1,j+1) - (U(i-1,j-1) + 2*U(i-1,j) + F2(i-1,j+1); Gy = (U(i-1,j+1) + 2*U(i,j+1) + F2(i+1,j+1) - (U(i-1,j-1) + 2*U(i,j-1) + F2(i+1,j-1); uSobel(i,j) = sqrt(Gx2 + Gy2); end end subplot(1,2,2);imshow(im2uint8(uSobel):title(Sobel邊緣檢

16、測后); （4） Laplace邊緣檢測主函數：f=imread(C:Users34899Desktoptimg.jpg)fGray=rgb2gray(f);figure()imshow(fGray),title(灰度圖像);fGray=double(fGray);T=60;LapModType=8;%設置laplace模板方式fLapEdge=LaplaceEdge(fGray,LapModType,T);fGrayLapEdge=uint8(fLapEdge);figure()imshow(fLapEdge),title(laplace邊緣圖像);四鄰域算子和八鄰域算子function f

17、LapEdge=LaplaceEdge(fGray,MoldType,Thresh)%-參數介紹-%輸入參數：% fGray;輸入的灰度圖像% MoldType:模板類型，包括四鄰域和八鄰域% Thresh:邊緣檢測閾值%輸出參數：% fEdge:邊緣像素點，存儲的是二值化圖像r,c=size(fGray);fLapEdge=zeros(r,c); %四鄰域拉普拉斯邊緣檢測算子if 4=MoldType for i=2:r-1 for j=2:c-1 Temp=-4*fGray(i,j)+fGray(i-1,j)+fGray(i+1,j)+fGray(i,j-1)+fGray(i,j+1);

18、if TempThresh fLapEdge(i,j)=255; else fLapEdge(i,j)=0; end end endend%八鄰域拉普拉斯邊緣檢測算子if 8=MoldType for i=2:r-1 for j=2:c-1 Temp=-8*fGray(i,j)+fGray(i-1,j)+fGray(i+1,j)+fGray(i,j-1)+fGray(i,j+1)+fGray(i-1,j-1)+fGray(i+1,j+1)+fGray(i+1,j-1)+fGray(i-1,j+1); if TempThresh fLapEdge(i,j)=255; else fLapEdge(

19、i,j)=0; end end endEnd（5） 拉普拉斯高斯檢測器clc; clear all; close all;hsize = 5 5;sigma = 0.5;h = fspecial(log, hsize, sigma);I = imread(C:Users34899Desktoptimg.jpg);bw = imfilter(I, h, replicate);figure;subplot(1, 2, 1); imshow(I, ); title(原圖像, FontWeight, Bold);subplot(1, 2, 2); imshow(bw, ); title(邊緣圖像, F

20、ontWeight, Bold);（6） Canny算子tic % 記錄CPU使用時間，開始計時I = imread(C:Users34899Desktoptimg.jpg); % 讀入圖像I=rgb2gray(I); % 轉化為灰色圖像% figure,imshow(I);title(原圖) % 顯示原圖像% a = im2single(I); % 將圖像矩陣轉換為單精度類型% a = im2uint8(a); % 將圖像矩陣轉換為8位無符號整數類型% I = uint16(I); % 8位無符號整數轉化為16位無符號整數 %* matlab的高斯平滑 *row,line = size(I)

21、; % 行列temp = zeros(row+8,line+8); % 申請空間imagedata = zeros(row,line); % 申請空間% 高斯核計算,核大小為9*9GaussianDieOff = .0001;sigma = 1;pw = 1:30; % possible widthsssq = sigma2;width = find(exp(-(pw.*pw)/(2*ssq)GaussianDieOff,1,last);if isempty(width) width = 1; % the user entered a really small sigmaendt = (-wi

22、dth:width);gau = exp(-(t.*t)/(2*ssq)/(2*pi*ssq); kernel = gau * gau; % 高斯核,大小為9*9kernel = im2single(kernel); % 圖像邊緣復制for q = 5:row+4 % 首先載入源圖像信息 for p = 5:line+4 temp(q,p) = I(q-4,p-4); endendfor q = 1:row+8 % 復制列5:line+4 for p = 5:line+4 if(qrow+4) temp(q,p) = temp(row+4,p); end endendfor q = 1:row

23、+8 % 復制行1:row+8 for p = 1:line+8 if(pline+4) temp(q,p) = temp(q,line+4); end endendfor q = 1:row % 高斯卷積 for p = 1:line for n = 1:9 for m = 1:9 imagedata(q,p) = imagedata(q,p) + kernel(n,m) * temp(q+n-1,p+m-1); end end endend%* matlab的高斯平滑說明結束 * %* matlab求梯度 *dx = zeros(size(I); % 申請空間,存放x方向的幅值dy = z

24、eros(size(I); % 申請空間,存放y方向的幅值x,y=meshgrid(-width:width,-width:width);dgau2D=-x.*exp(-(x.*x+y.*y)/(2*ssq)/(pi*ssq); % 核% 圖像邊緣復制for q = 5:row+4 % 首先載入源圖像信息 for p = 5:line+4 temp(q,p) = imagedata(q-4,p-4); endendfor q = 1:row+8 % 復制列5:line+4 for p = 5:line+4 if(qrow+4) temp(q,p) = temp(row+4,p); end en

25、dendfor q = 1:row+8 % 復制行1:row+8 for p = 1:line+8 if(pline+4) temp(q,p) = temp(q,line+4); end endendfor q = 1:row % x方向卷積 for p = 1:line for n = 1:9 for m = 1:9 dx(q,p) = dx(q,p) + dgau2D(n,m) * temp(q+n-1,p+m-1); end end endenddgau2D = dgau2D;for q = 1:row % y方向卷積 for p = 1:line for n = 1:9 for m =

26、 1:9 dy(q,p) = dy(q,p) + dgau2D(n,m) * temp(q+n-1,p+m-1); end end endend%* matlab求梯度結束 * %* 自適應閾值計算 *value = sqrt(dx.2 + dy.2); % 幅值value = double(value);value = value/max(max(value);% 歸一化 temp_value = value;counts=imhist(value, 64);high_threshold = find(cumsum(counts) 0.7*row*line,1,first) / 64;low

27、_threshold = 0.4 * high_threshold; % 低閾值clear counts%* 自適應閾值計算結束 * %* 方向與閾值判別 *num = 0; % 當前堆棧個數flag = zeros(80000,2); % 堆棧temp_flag = zeros(80000,2); % 臨時堆棧imagedata = zeros(row,line); % 初始化% 方向：% 0為045% 1為4590% 2為90135% 3為135180direction = zeros(size(I); % 申請空間,存放小于高閾值而大于低閾值for q = 2:row-1 for p =

28、 2:line-1 if(dy(q,p)-dy(q,p) | (dy(q,p)=0 & dx(q,p)0 & -dy(q,p)=dx(q,p) | (dx(q,p)0 & -dy(q,p)=dx(q,p) d = abs(dx(q,p)/dy(q,p); % 4590方向 gradmag = temp_value(q,p); gradmag1 = temp_value(q-1,p)*(1-d) + temp_value(q-1,p+1)*d; gradmag2 = temp_value(q+1,p)*(1-d) + temp_value(q+1,p-1)*d; elseif(dx(q,p)dy

29、(q,p) | (dx(q,p)=0 & dx(q,p)dy(q,p) d = abs(dx(q,p)/dy(q,p); % 90135方向 gradmag = temp_value(q,p); gradmag1 = temp_value(q-1,p)*(1-d) + temp_value(q-1,p-1)*d; gradmag2 = temp_value(q+1,p)*(1-d) + temp_value(q+1,p+1)*d; elseif(dy(q,p)0 & dx(q,p)0 & dx(q,p)=dy(q,p) d = abs(dy(q,p)/dx(q,p); % 135180方向 g

30、radmag = temp_value(q,p); gradmag1 = temp_value(q,p-1)*(1-d) + temp_value(q-1,p-1)*d; gradmag2 = temp_value(q,p+1)*(1-d) + temp_value(q+1,p+1)*d; end if(gradmag=gradmag1 & gradmag=gradmag2) if(gradmag = high_threshold) % 高閾值判別 value(q,p) = 255; elseif(gradmag = low_threshold)% 低閾值判別 value(q,p) = 125

31、; else value(q,p) = 0; % 小于低閾值 end else value(q,p) = 0; % 非極大值抑制 end end endfor q = 2:row-1 % 高閾值判別，檢查高閾值鄰域8個方向范圍內是否存在低閾值 for p = 2:line-1 if(value(q,p) = 255) imagedata(q,p) = 255; if(value(q-1,p-1) = 125) value(q-1,p-1) = 255; imagedata(q-1,p-1) = 255; if(q-1 1) & (p-1 1) num = num + 1; flag(num,1

32、) = q-1; flag(num,2) = p-1; end end if(value(q-1,p) = 125) value(q-1,p) = 255; imagedata(q-1,p) = 255; if(q-1 1) num = num + 1; flag(num,1) = q-1; flag(num,2) = p; end end if(value(q-1,p+1) = 125) value(q-1,p+1) = 255; imagedata(q-1,p+1) = 255; if(q-1 1) & (p+1 1) num = num + 1; flag(num,1) = q; fla

33、g(num,2) = p-1; end end if(value(q,p+1) = 125) value(q,p+1) = 255; imagedata(q,p+1) = 255; if(p+1 line) num = num + 1; flag(num,1) = q; flag(num,2) = p+1; end end if(value(q+1,p-1) = 125) value(q+1,p-1) = 255; imagedata(q+1,p-1) = 255; if(q+1 1) num = num + 1; flag(num,1) = q+1; flag(num,2) = p-1; e

34、nd end if(value(q+1,p) = 125) value(q+1,p) = 255; imagedata(q+1,p) = 255; if(q+1 row) num = num + 1; flag(num,1) = q+1; flag(num,2) = p; end end if(value(q+1,p+1) = 125) value(q+1,p+1) = 255; imagedata(q+1,p+1) = 255; if(q+1 row) & (p+1 1) & (p-1 1) num = num + 1; temp_flag(num,1) = q-1; temp_flag(n

35、um,2) = p-1; end end if(value(q-1,p) = 125) value(q-1,p) = 255; imagedata(q-1,p) = 255; if(q-1 1) num = num + 1; temp_flag(num,1) = q-1; temp_flag(num,2) = p; end end if(value(q-1,p+1) = 125) value(q-1,p+1) = 255; imagedata(q-1,p+1) = 255; if(q-1 1) & (p+1 1) num = num + 1; temp_flag(num,1) = q; temp_flag(num,2) = p-1; end end if(value(q,p+1) = 125) value(q,p+1) = 255; imagedata(q,p+1) = 255; if(p+1 line) num = num + 1; temp_flag(num,1) = q; temp_flag(num,2) = p+1;