%将im转换成灰度图和double型
im=rgb2gray(im);
im=im2double(im);
imshow(im),title('gray');
%由于噪点过多,用高斯模糊减少噪点
hsize=31;
sigma=10;
h=fspecial('gaussian',hsize,sigma);
img=imfilter(im,h);
%使用canny找到边界
img=edge(img,'canny');
figure,imshow(img),title('edge');
%求出图像大小
[row col]=size(img);
%设置每个格子的长宽
one_theta=1;
one_rho=1;
%计数矩阵
gridnum_theta=ceil(180/one_theta); %theta轴的格子数180,坐标点的数量由于有0所以要再+1
gridnum_rho=(ceil(sqrt(row^2+col^2)/one_rho))*2+1; % ceil(sqrt(row^2+col^2)/one_rho)是rho正半轴的值,所以整个rho点的个数为其*2+1
accum=zeros(gridnum_rho,gridnum_theta+1); %得到投票用的空矩阵
p=ceil(sqrt(row^2+col^2)/one_rho); %p+1处为坐标为0处
%投票
for i=1:row %遍历每个img上的像素点
for j=1:col
if img(i,j)==1 %如果该点是边界点,则遍历θ
for theta_s=-90:90 %遍历θ,计算rho,在每个(rho,theta)的点加一票
rho_s=ceil(j*cosd(theta_s)+i*sind(theta_s));
accum((p+1+rho_s),theta_s+90+1)=accum((p+1+rho_s),theta_s+90+1)+1; %这里要注意accum的坐标索引要为正数,坐标应要包含原点
end
end
end
end
%建立accum的rho和theta两个坐标轴矩阵(都是1*n的)
rho=zeros(1,length(gridnum_rho));
for i=0:2*p+1
rho(1,i+1)=p*(-1)+i;
end
theta=zeros(1,length(gridnum_theta)+1);
for i=0:180+1
theta(1,i+1)=-90+i;
end
%找到最大的n个点
n=30;
points=zeros(n,2);
accum_temp=accum; %复制一个accum,然后在其中找最大值,每找到一个记录下坐标,就将该值改为0。
k=1;
for i2=1:n
m=max(max(accum_temp(:)));
for i=1:gridnum_rho
for j=1:gridnum_theta+1
if accum_temp(i,j)==m
if m==0
return
end
points(k,1)=i;
points(k,2)=j;
accum_temp(i,j)=0;
k=k+1;
if k>n
return
end
end
end
end
end
points_rhotheta=zeros(n,2);
%找到对应极坐标位置
for i=1:n
if points(i,1)==0 & points(i,2)==0 %有可能没有取满peaks的n个点,记录peaks点坐标的矩阵会有(0,0)
return
end
points_rhotheta(i,1)=rho(1,points(i,1));
points_rhotheta(i,2)=theta(1,points(i,2));
end
%画线
figure,imshow(im),title('line');
hold on;
x=linspace(1,col);
y=linspace(1,row);
for i=1:n
y=(points_rhotheta(i,1)-x*cosd(points_rhotheta(i,2)))/sind(points_rhotheta(i,2));
plot(x,y,'LineWidth',2,'Color','green');
end
hold off;