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

2、分構(gòu)成：2fx,y=2f(x,y)x2+2f(x,y)y22. 邊緣檢測(cè)的基本思想：（1） 尋找灰度的一階導(dǎo)數(shù)的幅度大于某個(gè)指定閾值的位置；（2） 尋找灰度的二階導(dǎo)數(shù)有零交叉的位置。3. 幾種方法簡(jiǎn)介（1） Sobel邊緣檢測(cè)器：以差分來(lái)代替一階導(dǎo)數(shù)。Sobel邊緣檢測(cè)器使用一個(gè)33鄰域的行和列之間的離散差來(lái)計(jì)算梯度，其中，每行或每列的中心像素用2來(lái)加權(quán)，以提供平滑效果。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邊緣檢測(cè)器:使用下圖所示模板來(lái)數(shù)字化

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

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

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

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

7、算子的形式，用于提取邊緣，可以利用快速卷積函數(shù)， 簡(jiǎn)單有效，因此應(yīng)用廣泛。美中不足的是，Sobel算子并沒(méi)有將圖像的主體與背景嚴(yán)格地區(qū)分開(kāi)來(lái)，即Sobel算子沒(méi)有嚴(yán)格地模擬人的視覺(jué)生理特征，所以提取的圖像輪廓有時(shí)并不能令人滿(mǎn)意。（4） Laplace算子：Laplace算子是一種各向同性算子，在只關(guān)心邊緣的位置而不考慮其周?chē)南笏鼗叶炔钪禃r(shí)比較合適。Laplace算子對(duì)孤立象素的響應(yīng)要比對(duì)邊緣或線(xiàn)的響應(yīng)要更強(qiáng)烈，因此只適用于無(wú)噪聲圖象。存在噪聲情況下，使用Laplacian算子檢測(cè)邊緣之前需要先進(jìn)行低通濾波。（5） LoG算子：由于二階導(dǎo)數(shù)是線(xiàn)性運(yùn)算，利用LoG卷積一幅圖像與首先使用高斯型平滑

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

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

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

12、ourcePic=imread(C:Users34899Desktoptimg.jpg); %讀取原圖像grayPic=mat2gray(sourcePic);%實(shí)現(xiàn)圖像矩陣的歸一化操作m,n=size(grayPic);newGrayPic=grayPic;%為保留圖像的邊緣一個(gè)像素PrewittNum=0;%經(jīng)Prewitt算子計(jì)算得到的每個(gè)像素的值PrewittThreshold=0.5;%設(shè)定閾值for j=2:m-1 %進(jìn)行邊界提取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邊緣檢測(cè) clear all; close all; imag = imread(C:Users34899Desktoptimg.jpg); %讀取關(guān)鍵幀 imag = rgb2gray(imag); %轉(zhuǎn)化為灰度圖 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邊緣檢測(cè) 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、測(cè)后); （4） Laplace邊緣檢測(cè)主函數(shù)：f=imread(C:Users34899Desktoptimg.jpg)fGray=rgb2gray(f);figure()imshow(fGray),title(灰度圖像);fGray=double(fGray);T=60;LapModType=8;%設(shè)置laplace模板方式fLapEdge=LaplaceEdge(fGray,LapModType,T);fGrayLapEdge=uint8(fLapEdge);figure()imshow(fLapEdge),title(laplace邊緣圖像);四鄰域算子和八鄰域算子function f

17、LapEdge=LaplaceEdge(fGray,MoldType,Thresh)%-參數(shù)介紹-%輸入?yún)?shù)：% fGray;輸入的灰度圖像% MoldType:模板類(lèi)型，包括四鄰域和八鄰域% Thresh:邊緣檢測(cè)閾值%輸出參數(shù)：% fEdge:邊緣像素點(diǎn)，存儲(chǔ)的是二值化圖像r,c=size(fGray);fLapEdge=zeros(r,c); %四鄰域拉普拉斯邊緣檢測(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%八鄰域拉普拉斯邊緣檢測(cè)算子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） 拉普拉斯高斯檢測(cè)器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使用時(shí)間，開(kāi)始計(jì)時(shí)I = imread(C:Users34899Desktoptimg.jpg); % 讀入圖像I=rgb2gray(I); % 轉(zhuǎn)化為灰色圖像% figure,imshow(I);title(原圖) % 顯示原圖像% a = im2single(I); % 將圖像矩陣轉(zhuǎn)換為單精度類(lèi)型% a = im2uint8(a); % 將圖像矩陣轉(zhuǎn)換為8位無(wú)符號(hào)整數(shù)類(lèi)型% I = uint16(I); % 8位無(wú)符號(hào)整數(shù)轉(zhuǎn)化為16位無(wú)符號(hào)整數(shù) %* matlab的高斯平滑 *row,line = size(I)

21、; % 行列temp = zeros(row+8,line+8); % 申請(qǐng)空間imagedata = zeros(row,line); % 申請(qǐng)空間% 高斯核計(jì)算,核大小為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); % 圖像邊緣復(fù)制for q = 5:row+4 % 首先載入源圖像信息 for p = 5:line+4 temp(q,p) = I(q-4,p-4); endendfor q = 1:row+8 % 復(fù)制列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 % 復(fù)制行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的高斯平滑說(shuō)明結(jié)束 * %* matlab求梯度 *dx = zeros(size(I); % 申請(qǐng)空間,存放x方向的幅值dy = z

24、eros(size(I); % 申請(qǐng)空間,存放y方向的幅值x,y=meshgrid(-width:width,-width:width);dgau2D=-x.*exp(-(x.*x+y.*y)/(2*ssq)/(pi*ssq); % 核% 圖像邊緣復(fù)制for q = 5:row+4 % 首先載入源圖像信息 for p = 5:line+4 temp(q,p) = imagedata(q-4,p-4); endendfor q = 1:row+8 % 復(fù)制列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 % 復(fù)制行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求梯度結(jié)束 * %* 自適應(yīng)閾值計(jì)算 *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%* 自適應(yīng)閾值計(jì)算結(jié)束 * %* 方向與閾值判別 *num = 0; % 當(dāng)前堆棧個(gè)數(shù)flag = zeros(80000,2); % 堆棧temp_flag = zeros(80000,2); % 臨時(shí)堆棧imagedata = zeros(row,line); % 初始化% 方向：% 0為045% 1為4590% 2為90135% 3為135180direction = zeros(size(I); % 申請(qǐng)空間,存放小于高閾值而大于低閾值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個(gè)方向范圍內(nèi)是否存在低閾值 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;