canny算子的matlab程序

function canny()
clc
clear all
I = imread('lingjian.JPG');%读入图像
I=imresize(I,[480,640]);
imshow(I);
title('原图');
I1 = rgb2gray(I);%变换成灰度图像
imshow(I1);
BW =edge(I1,'canny'); % 调用canny函数
figure,imshow(BW); % 显示分割后的图像，即梯度图像
title('Canny')
[n,m] = size(BW);
a_max=40; %a为椭圆长轴，范围需要自己根据图片中的椭圆形状确定
a_min=18;
b_max=30; %b为椭圆短轴
b_min=4;
step_a=3; %步长
step_b=3;
theta_min=-pi/18; %theta倾斜角
theta_max=pi/9;
step_theta=pi/60;
step_angle=2;
yz =0.82; %阈值（自己定）

size_a = round((a_max-a_min)/step_a)+1;
size_b = round((b_max-b_min)/step_b)+1;
size_theta=round((theta_max-theta_min)/step_theta)+1;
size_angle = round(2*pi/step_angle);

hough_space = zeros(m,n,size_a,size_b,size_theta); %设定5维的参数空间并使初值为0

[cols,rows] = find(BW); %find—找出非零元素的索引
ecount = size(rows);

% Hough
%图像空间(x,y)对应到参数空间(a,b,p、q、theta)
% p = x-a*cos(angle)*cos(theta)+b*sin(angle)*sin(theta)
% q = y-a*cos(angle)*sin(theta)-b*sin(angle)*cos(theta)
%rows(i)行坐标
for i=1:ecount for a=1:size_a
for b=1:size_b
for theta=1:size_theta
for k=1:size_angle
%hough变换
p = round(rows(i)-(a_min+(a-1)*step_a)*cos(k*step_angle)*cos(theta_min+(theta-1)*step_theta)+(b_min+(b-1)*step_b)*sin(k*step_angle)*sin(theta_min+(theta-1)*step_theta));
q = round(cols(i)-(a_min+(a-1)*step_a)*cos(k*step_angle)*sin(theta_min+(theta-1)*step_theta)-(b_min+(b-1)*step_b)*sin(k*step_angle)*cos(theta_min+(theta-1)*step_theta));
if(p>0&p<=m&q>0&q<=n)
hough_space(p,q,a,b,theta) = hough_space(p,q,a,b,theta)+1;
end
end
end
end
end
end

% 搜索超过阈值的聚焦点
max_para = max(max(max(max(max(hough_space)))));
index = find(hough_space>max_para*yz); %find—找出hough_space中大于阈值的缩引并存入 index
length = size(index);
hough_circle1=zeros(m,n); %确定为椭圆上的点的坐标
hough_circle2=zeros(m,n);

%找出峰值对应的参数空间坐标
for k=1:length
par5 = floor((index(k)-1)/(m*n*size_a*size_b))+1; %theta增量
par4 = floor((index(k)-(par5-1)*(m*n*size_a*size_b))/(m*n*size_a))+1; %b增量
par3 = floor((index(k)-(par5-1)*(m*n*size_a*size_b)-(par4-1)*(m*n*size_a))/(m*n))+1; %a 增量
par2 = floor((index(k)-(par5-1)*(m*n*size_a*size_b)-(par4-1)*(m*n*size_a)-(par3-1)*(m*n))/m)+1; %p增量
par1 = index(k)-(par5-1)*(m*n*size_a*size_b)-(par4-1)*(m*n*size_a)-(par3-1)*(m*n)-(par2-1)*m; %q增量
par5=theta_min+(par5-1)*step_theta;
par4 = b

_min+(par4-1)*step_b;
par3 = a_min+(par3-1)*step_a;
theta(k)=par5;
b(k)=par4;
a(k)=par3;
q(k)=par2;
p(k)=par1;
end

%求出两圆参数平均值
[row1 col1]=size(p);
count=1;
theta=sort(theta);
p=sort(p);
q=sort(q);
a=sort(a);
b=sort(b);
THETA(count)=theta(1);
P(count)=p(1);
A(count)=a(1);
B(count)=b(1);
Q(count)=q(1);
for t=1:1:col1
if abs(P(count)-p(t))<=10
THETA(count)=(theta(t)+THETA(count))/2;
A(count)=(a(t)+A(count))/2;
B(count)=(b(t)+B(count))/2;
P(count)=(p(t)+P(count))/2;
Q(count)=(q(t)+Q(count))/2;
else
count=count+1;
THETA(count)=theta(t);
A(count)=a(t);
B(count)=b(t);
P(count)=p(t);
Q(count)=q(t);
end
end
THETA
A
B
P
Q
%绘制椭圆
TYZ=zeros(1,2);
TYY=zeros(1,2);
ct_z=1;
ct_y=1;
for i=1:ecount
if round(((rows(i)-P(1))*cos(THETA(1))+(cols(i)-Q(1))*sin(THETA(1)))^2/(A(1)^2)+(-(rows(i)-P(1))*sin(THETA(1))+(cols(i)-Q(1))*cos(THETA(1)))^2/(B(1)^2))<1.5 ...
&round(((rows(i)-P(1))*cos(THETA(1))+(cols(i)-Q(1))*sin(THETA(1)))^2/(A(1)^2)+(-(rows(i)-P(1))*sin(THETA(1))+(cols(i)-Q(1))*cos(THETA(1)))^2/(B(1)^2))>0.5
TYY(ct_y,1)=rows(i);
TYY(ct_y,2)=cols(i);
hough_circle1(cols(i),rows(i))=1;
ct_y=ct_y+1;
end
if round(((rows(i)-P(2))*cos(THETA(2))+(cols(i)-Q(2))*sin(THETA(2)))^2/(A(2)^2)+(-(rows(i)-P(2))*sin(THETA(2))+(cols(i)-Q(2))*cos(THETA(2)))^2/(B(2)^2))<1.5 ...
&round(((rows(i)-P(2))*cos(THETA(2))+(cols(i)-Q(2))*sin(THETA(2)))^2/(A(2)^2)+(-(rows(i)-P(2))*sin(THETA(2))+(cols(i)-Q(2))*cos(THETA(2)))^2/(B(2)^2))>0.5
TYZ(ct_z,1)=rows(i);
TYZ(ct_z,2)=cols(i);
hough_circle2(cols(i),rows(i))=1;
ct_z=ct_z+1;
end
end
figure
imshow(hough_circle1),title('hough圆1')
figure
imshow(hough_circle2),title('hough圆2')

%分别计算两个椭圆的参数
[row_TYZ,col_TYZ]=size(TYZ);
for i1=1:1:row_TYZ
WTZ(i1,:)=[TYZ(i1,1)^2 TYZ(i1,1)*TYZ(i1,2) TYZ(i1,2)^2 TYZ(i1,1) TYZ(i1,2) 1; ];
end
[v_z1,d_z1]=svd(WTZ'*WTZ);
v_z1=vpa(v_z1,8)
d_z1=double(d_z1);

[row_TYY,col_TYY]=size(TYY);
for j1=1:1:row_TYY
WTY(j1,:)=[TYY(j1,1)^2 TYY(j1,1)*TYY(j1,2) TYY(j1,2)^2 TYY(j1,1) TYY(j1,2) 1; ];
end
[v_y1,d_y1]=svd(WTY'*WTY);
v_y1=vpa(v_y1,8)
d_y1=double(d_y1);