clc;
clear all;
close all;
%边缘连接测试图像
I=im2double(imread('D:\Gray Files\10-34.tif'));
[M,N]=size(I);
n=13;
%寻找图像的边缘
sigma=2;
H=0.15;
L=0.05;
g_c=CannyEdgeDetector(n,sigma,H,L,I);
%找出所有的边线点
[rows,cols]=find(g_c);
points=cat(2,rows,cols);
%进行霍夫变换，Hough Transform
theta_max=90;
theta_range=-theta_max:theta_max;
n=length(theta_range);
rho_max=floor(sqrt(M^2+N^2));
rho_range=-rho_max:rho_max;
m=length(rho_range);
Hough=zeros(m,n);
%霍夫变换后的矩阵,行与上述边界点集相对应，列为每一个点的霍夫转换
while ~isempty(points)
%取出堆栈顶上的点
p=points(1,:);
%霍夫变换
for k=-90:1:90
j=k+theta_max+1;
%四舍五入
i=round(cosd(k)*p(1,1)+sind(k)*p(1,2))+rho_max+1;
Hough(i,j)=Hough(i,j)+1;
end
%删除该点
points(1,:)=[];
end
Hou=Hough/max(Hough(:));
for i=1:m
for j=1:n
if Hou(i,j)>0
Hou(i,j)=Hou(i,j)+0.5;
end
end
end
%显示边界线的霍夫空间图像
% imshow(Hou)
%寻找最大值及其在霍夫空间的坐标
%霍夫空间门限
th=200;
T=max(Hough(:));
[rows,cols]=find(Hough==T)
inds=cat(2,rows,cols);
%求解笛卡尔坐标中的线
ps=[];
while ~isempty(inds)
%取出极坐标值
ind=inds(1,:);
rho=ind(1,1);
theta=ind(1,2);
%还原角度和极坐标值
rho=rho-rho_max-1;
theta=theta-theta_max-1;
%定义一条线
g_l=zeros(M,N);
%判断theta的角度
if theta==0
%取出边缘图像中第rho行为1的点
ind=find(g_c(rho,:)==1);
%取出最小下标和最大下标
ind_y_min=ind(1);
ind_y_max=ind(end);
g_l(rho,ind_y_min:ind_y_max)=1;
%与边缘图像交集
g_c=g_c+g_l;
else
for i=1:M
j=round((rho-i*cosd(theta))/sind(theta));
if j>0 && j<=N
g_l(i,j)=1;
end
end
g_l_c=and(g_l,g_c);
[rows,cols]=find(g_l_c==1);
ps=cat(2,rows,cols);
ps=sortrows(ps,1);
rows=sort(rows);
ind_x_min=rows(1,1);
ind_x_max=rows(end);
if ind_x_min>1
g_l(1:ind_x_min-1,:)=0;
end
if ind_x_max<M
g_l(ind_x_max+1:M,:)=0;
end
g_c=g_c+g_l;
end
inds(1,:)=[];
end
imshow(g_c)
hold on
plot(ps(:,2),ps(:,1),'red','LineWidth',2)
hold off
Canny算子，CannyEdgeDetector如下：
%canny边界检测器
% n 高斯低通滤波器的大小
% sigma 高斯滤波器参数
% H 高门限
% L 低门限
function [g]=CannyEdgeDetector(n,sigma,H,L,I)
[M,N]=size(I);
%%
%=============================边缘检测(五)=================================
% Canny Edge Detector
%-------------------------用高斯低通滤波器平滑图像-------------------------
%建筑图像所设参数
% n=25;
% sigma=4;
%头颅扫描图像所设参数
% n=13;
% sigma=2;
%lena测试图像所设参数
% n=5;
% sigma=1;
type='symmetric';
f_s=GaussianBlur(n,I,sigma,type);
% imshow(f_s)
% imwrite(f_s,'D:\Gray Files\lena-test.jpg','jpg');
%-----------------------计算平滑后的图像梯度和角度-------------------------
n_l=1;
%Sobel算子
s_y=[-1 -2 -1;
0 0 0;
1 2 1];
s_x=[-1 0 1;
-2 0 2;
-1 0 1];
%定义梯度和角度
gx=zeros(M,N);
gy=zeros(M,N);
f_s_pad=padarray(f_s,[n_l,n_l],'replicate');
for i=1:M
for j=1:N
Block=f_s_pad(i:i+2*n_l,j:j+2*n_l);
gx(i,j)=sum(sum(Block.*s_x));
gy(i,j)=sum(sum(Block.*s_y));
end
end
type='replicate';
gx=GaussianBlur(n,gx,sigma,type);
gy=GaussianBlur(n,gy,sigma,type);
M_s=sqrt(gx.^2+gy.^2);
M_s=M_s/max(M_s(:));
a_s=atan2(gy,gx)*180/pi;
%--------------------对梯度图像进行差值非极大值抑制------------------------
n_l=1;
%定义非极大值抑制图像
g_N=M_s;
M_s_pad=padarray(M_s,[n_l,n_l],'replicate');
for i=1:M
for j=1:N
%取出中心点的梯度值
K=M_s_pad(i+1,j+1);
theta=a_s(i,j);
if (theta>=0 && theta<=45) ||...
(theta<-135 && theta>=-180)
yBot=[M_s_pad(i+1,j+2) M_s_pad(i+2,j+2)];
yTop=[M_s_pad(i+1,j) M_s_pad(i,j)];
k=abs(gy(i,j)/M_s_pad(i+1,j+1));
K1=(yBot(2)-yBot(1))*k+yBot(1);
K2=(yTop(2)-yTop(1))*k+yTop(1);
end
if (theta>45 && theta<=90) ||...
(theta<-90 && theta>=-135)
yBot=[M_s_pad(i+2,j+1) M_s_pad(i+2,j+2)];
yTop=[M_s_pad(i,j+1) M_s_pad(i,j)];
k=abs(gx(i,j)/M_s_pad(i+1,j+1));
K1=(yBot(2)-yBot(1))*k+yBot(1);
K2=(yTop(2)-yTop(1))*k+yTop(1);
end
if (theta>90 && theta<=135) ||...
(theta<-45 && theta>=-90)
yBot=[M_s_pad(i+2,j+1) M_s_pad(i+2,j)];
yTop=[M_s_pad(i,j+1) M_s_pad(i,j+2)];
k=abs(gx(i,j)/M_s_pad(i+1,j+1));
K1=(yBot(2)-yBot(1))*k+yBot(1);
K2=(yTop(2)-yTop(1))*k+yTop(1);
end
if (theta>135 && theta<=180) ||...
(theta<0&& theta>=-45)
yBot=[M_s_pad(i+1,j) M_s_pad(i+2,j)];
yTop=[M_s_pad(i+1,j+2) M_s_pad(i,j+2)];
k=abs(gy(i,j)/M_s_pad(i+1,j+1));
K1=(yBot(2)-yBot(1))*k+yBot(1);
K2=(yTop(2)-yTop(1))*k+yTop(1);
end
if K<K1 || K<K2
g_N(i,j)=0;
end
end
end
g_N=g_N/max(g_N(:));
imshow(g_N)
%-------------------------------设置双门限---------------------------------
%对非极大值抑制图像进行扩展(0)，方便后续的边界处理
n_l=1;
g_N_pad=padarray(g_N,[n_l,n_l],'replicate');
T_max=max(g_N_pad(:));
%lena图像所设参数
% T_max=max(g_N_pad(:));
% H=0.275;
% L=0.25;
%高门限
% T_H=T_max*H;
%低门限
% T_L=T_H*L;
%头颅图像所设参数
% T_H=0.15;
% T_L=0.05;
T_H=H;
T_L=L;
%建筑图像所设参数
% H=0.1;
% L=0.04;
%高门限
% T_H=T_max*H;
% %低门限
% T_L=T_H*L;
%寻找大于高门限点
ind_H=find(g_N_pad>T_H);
g_N_pad(ind_H)=1;
ind_Zeros=find(g_N_pad<T_L);
g_N_pad(ind_Zeros)=0;
% imshow(g_N_pad)
%--------------------------消除未连通的边界点------------------------------
M_temp=M+2;
g_N_copy=g_N_pad;
g_N_copy(ind_H)=0;
ind_L=find(g_N_copy~=0);
g_N_pad(ind_L)=1;
while ~isempty(ind_H)
%取出第一个点
p=ind_H(1,1);
%初始化连通点集
Conn=[];
%查找p点是否有连通点，即在ind_L中是否有值
[ind_L,Conn]=FindConnected_8(p,ind_L,Conn,M_temp);
%向下查找所有与p点连通的点
while ~isempty(Conn)
%从连通集中取出第一个点，而后删除该点
p1=Conn(1,:);
Conn(1,:)=[];
[ind_L,Conn]=FindConnected_8(p1,ind_L,Conn,M_temp);
end
%标记p点已访问，即从ind_H中删除
ind_H(1,:)=[];
end
%将ind_L中未与ind_H连通的点，在图像中置零
g_N_pad(ind_L)=0;
g_N=g_N_pad(2:M+1,2:N+1);
%-----------------------------图像细化处理-----------------------------
[g]=ImageThinning(g_N);
% imshow(g)
end
Canny算子中的其他函数，如GaussianBlur、FindConnected_8、ImageThinning参见上一篇博文，这里就不在重复粘贴了