边缘提取，大津算法
?
?
解：采用Canny邊緣檢測器，代碼如下：
?
function?Canny()%自己寫的Canny算法，陳焜
clc;
clear?all;
close?all;
img=imread('Fig1006(a)(building).tif');
figure(1);subplot(221);imshow(img);title('原圖像');
figure(11);imshow(img);title('原圖像');
?
img=double(img);
[m,n]=size(img);
?
%%高斯低通濾波，起平滑作用
img1=GLPF(img);%自己寫的高斯低通濾波器程序?
figure(1);subplot(222);imshow(uint8(img1));title('高斯濾波');
figure(12);imshow(uint8(img1));title('高斯濾波后的圖像');
?
%%計算梯度幅值
Sobelx=[-1,-2,-1;0,0,0;1,2,1];%x方向的Sobel模板
Sobely=[-1,0,1;-2,0,2;-1,0,1];%y方向的Sobel模板
gx=imfilter(img1,Sobelx,'replicate');?%求x方向梯度分量,replicate表示圖像大小通過復制外邊界來擴展
gy=imfilter(img1,Sobely,'replicate');?%求y方向梯度分量
Mg=sqrt(gx.^2+gy.^2);?%梯度幅值
Ag=atan(gy./gx);?????%梯度方向
?
?
%%非最大抑制
img2=ImaxInhibit(Mg,Ag);
figure(15);imshow(uint8(img2));title('非最大抑制');
figure(1);subplot(223);imshow(uint8(img2));title('非最大抑制');
?
?
%%雙閾值（滯后閾值）處理
img3=DualThreshold(img2);
figure(20);imshow(img3);title('Canny處理圖像');
figure(1);subplot(224);imshow(img3);title('Canny處理圖像');
end
?
?
?
%----------------------------------------------------------------------
%以下為子函數
?
%高斯低通濾波子函數
function?gauss=GLPF(I)?
[M,?N]=size(I);
P=2*M;%填充圖像為P*Q
Q=2*N;
F0=zeros(P,Q);
F0(1:M,1:N)=I;?%填充圖像
F1=?fftshift(F0);%將頻域原點移到圖像中心；
F=?fft2(F1);?%?傅立葉變換
D0?=?400;?%?D0為截止頻率
for?u=1:P
????for?v=1:Q
????????D(u,v)=?sqrt((u-P/2)^2+(v-Q/2)^2);?%距離
????????H(u,v)=?1-exp(-D(u,v)^2/(2*D0^2));??%?Gauss濾波器函數
%????????G(u,v)?=?H(u,v).*F(u,v);
????end
end
G?=?H.*F;%放在循環外，減少計算量
g?=?ifft2(G);%?傅立葉反變換
g1=?real(g);?%取實部
g2=?ifftshift(g1);%將頻域原點移到圖像中心；
gauss=g2(1:M,1:N);%裁截圖像
end
?
?
?
?
?
?
?
?
?
%%非最大抑制子函數
function?Img=ImaxInhibit(Mg,Ag)
[M,N]=size(Mg);
Ag=Ag*180;
Img=Mg;
for?i=2:M-1
????for?j=2:N-1
???????if?Ag(i,j)>=-22.5&&Ag(i,j)<22.5||abs(Ag(i,j))>=157.5?%水平邊緣
?????????????if?Mg(i,j)<Mg(i+1,j)||Mg(i,j)<Mg(i-1,j)
????????????????Img(i,j)=0;?%抑制
?????????????end
????????Else
If
Ag(i,j)>=-67.5&&Ag(i,j)<-22.5||Ag(i,j)>=112.5&&Ag(i,j)<157.5%+45°
????????????????if?Mg(i,j)<Mg(i+1,j-1)||Mg(i,j)<Mg(i-1,j+1)
????????????????????Img(i,j)=0;?%抑制
????????????????end
????????else?
if?
Ag(i,j)>=-112.5&&Ag(i,j)<-67.5||Ag(i,j)>=67.5&&Ag(i,j)<112.5%垂直
????????????????if?Mg(i,j)<Mg(i,j-1)||Mg(i,j)<Mg(i,j+1)
????????????????????????Img(i,j)=0;?%抑制
????????????????end
?????????else?
if?
Ag(i,j)>-157.5&&Ag(i,j)<-112.5||Ag(i,j)>=22.5&&Ag(i,j)<67.5?%—45°
????????????????if?Mg(i,j)<Mg(i-1,j-1)||Mg(i,j)<Mg(i+1,j+1)
????????????????????????????Img(i,j)=0;?%抑制
????????????????end
??????????end
??????????end
??????????end
??????????end
????end
?end
end
?
?
?
?
?
?
?
?
%%雙閾值（滯后閾值）處理子函數
function?G=DualThreshold(Gn)
[M,N]=size(Gn);
Gn=mat2gray(Gn);%歸一化
?
Gnh=Gn;
Gnl=Gn;
Th=0.12;
Tl=0.04;%滿足高閾值與低閾值之比為3:1或者2:1
for?i=2:M-1
????for?j=2:N-1
????????if?Gn(i,j)<Th
????????????Gnh(i,j)=0;%高閾值抑制
????????end
????????if?Gn(i,j)<Tl
????????????Gnl(i,j)=0;%低閾值抑制
????????end
????end
end
?
Gnl1=Gnl-Gnh;%從Gnl中刪除所有來自Gnh的非零像素，參見教材P465?式（10.2-35）
?
B=[0,1,0;1,1,1;0,1,0];%用膨脹的方法進行連接
Gnl1=imdilate(Gnl1,B);
?
G=Gnh+Gnl1;%將Gnl1中的所有非零像素附加到Gnh中
G=Gnh;
?
figure(18);imshow(Gnl);title('Gnl');
figure(19);imshow(Gnh);title('Gnh');
figure(17);imshow(Gnl1);title('Gnl-Gnh');
end
?
?
?
?
?
?
?
?
?
?
?
?
運行結果如下:
?
?
?
?
?
?
?
?
?
?
解：大津算法，代碼如下：
?
function?Otsu()%自己寫的大津算法，ck
clc;
clear?all;
close?all;
Img=imread('Fig1013(a)(scanned-text-grayscale).tif');?%讀入圖片
Img=rgb2gray(Img);%因為原圖是3維的，故需要轉換
figure(1),imshow(Img);title('原圖像')
?
p=zeros(1,256);%存放各灰度級的比率
Img=double(Img);%雙精度化
[M,N]=size(Img);
%計算圖像直方圖???
for?k=0:255
????p(k+1)=length(find(Img==k))/(M*N);?%計算每級灰度出現的概率(數組下標從1開始)
end???????????
k=1:1:256;
figure(2);?bar(k,p);?title('原圖像直方圖');
?
for?i=2:256
????????if?p(i)~=0
????????????minT=i+1;%尋找比率不為0的最小灰度值
????????????break
????????end
end
?
for?i=256:-1:1
????????if?p(i)~=0;
????????????maxT=i-1;%找出比率不為0的最大灰度值
????????????break
????????end
end
?
Mg=0;%Mg為全局均值
for?i=1:255
??Mg=Mg+i*p(i);
end
?
Var=zeros(1,(maxT-minT));%類間方差
Var1=0;%最大類間方差
for?k=minT:maxT
?????P1=0;%被分到C1的概率
?????m=0;%被分到C1的像素的平均灰度值
????for?i=1:k
????????P1=P1+p(i);%求被分到C1的概率，?參照教材P481?式（10.3-4）
????????m=m+i*p(i);%求被分到C1的像素的平均灰度值，?參照教材P481?式（10.3-8）???????
????end
????Var(k)=(Mg*P1-m)^2/(P1*(1-P1));%求類間方差，?參照教材P481?式（10.3-15）
????
????if(Var(k)>=Var1)
???????Var1=Var(k);%取得最大類間方差
????end
end
?
%如果有多個k使類間方差取得最大值，則取多個k的平均值作為最佳閾值，參照教材P481
K=zeros(1,(maxT-minT));%最佳閾值
j=0;%統計取得最大類間方差的閾值個數
for?k=minT:maxT
????if(Var(k)==Var1)
???????j=j+1;
???????K(j)=k;%取得最大類間方差的閾
????end
end
K1=sum(K)/j?%取多個k的平均值作為最佳閾值
?
%用Ostu方法的最佳閾值處理
for?i=1:M
????????for?j=1:N
????????????if?(Img(i,j)>K1)
???????????????Img1(i,j)=255;
????????????else
???????????????Img1(i,j)=0;
????????????end
????????end
End
?
figure(3);imhist(Img1);title('Otsu直方圖');
figure(4);imshow(Img1);title('Otsu最佳閾值處理圖像');
text('Position',[400,480],'String','Otsu最佳閾值:','color','r')%標注最佳閾值
text('Position',[510,480],'String',num2str(K1),'color','r')
?
運行結果如下：
?
?
?
?
?
?
?
?
?
?
最優全局閾值：102
?
?
?
?
轉載于:https://www.cnblogs.com/chenkun1/p/4153689.html
總結
- 上一篇： 【转帖】计算机世界：后DRM时代的数字音
- 下一篇： job 做 ha 问题?