边缘提取，大津算法
- 试编程实现提取图像Fig1006(a)(building).tif中的边缘.
解：采用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
运行结果如下:
- 试编程实现图像最优全局阈值分割方法（即global optimal segmentation，大津算法），并将其应用于分割出图像Fig1013(a)(scanned-text-grayscale).tif中的文字，并求出最优全局阈值（global optimal threshold）.
解：大津算法，代码如下：
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
