% houghtransform sample program by Gerox (c) 2004 
close all;
clear all
xmax=120;
ymax=120;
rhomax=120;
thetamax=180;
%%%%%%%%%Define line parameter%%%%%%%%%%%%%%%%%%
%for rho=-50:10:50
rho=20;
%for theta=(-60:10:60)*pi/180;
theta=60*pi/180;
%%%%%%%%%make image%%%%%%%%%%%%%%%%%%%
bw=zeros(ymax,xmax);
if(abs(theta-pi/2)<pi/4)
	x=[-xmax/2 xmax/2-1];
	y=rho/sin(theta)-x*cos(theta)/sin(theta);
	xindex=x+1+xmax/2;
	yindex=round(ymax/2-rho/sin(theta)+(xindex-1-xmax/2)*cos(theta)/sin(theta));
else
	y=[ymax/2-1 -ymax/2];
	x=rho/cos(theta)-y*sin(theta)/cos(theta);
	yindex=-y+ymax/2;
	xindex=round(rho/cos(theta)-(ymax/2-yindex)*sin(theta)/cos(theta)+xmax/2+1);
end
bw=matline(bw,yindex(1),xindex(1),yindex(2),xindex(2));
subplot(2,2,1);imshow(~bw);
title(['原画像 \rho=',num2str(rho),' \theta=',num2str(theta*180/pi)]);
%%%%%%%%Hough transform %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hbw=houghf(bw,rhomax,thetamax);
%%%%%%%%Calculate rho and theta from Hough pace %%%%%%%%%%%%%%%%%%%%%%%%%
[ymax xmax]=size(bw);% This is not necessary since these varibales are already defined.
[rhoindex,thetaindex]=find(hbw==max(hbw(:)));
rhoindex=mean(rhoindex);
thetaindex=mean(thetaindex);
rho=(rhomax/2-rhoindex)*sqrt(xmax*xmax+ymax*ymax)/rhomax
theta=pi/thetamax*(thetaindex-1)
[rho theta*180/pi]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
subplot(2,2,3);surf(((1:thetamax)/thetamax)*180,-((1:rhomax)/rhomax-0.5)*sqrt(xmax*xmax+ymax*ymax),log(hbw));axis tight;colormap hot;shading interp;view(2);axis square
title('Hough 空間でのイメージ');ylabel('\rho');xlabel('\theta');colormap hot
%%%%%%%Reconstruct line from calculated rho and theta %%%%%%%%%%%%%%%%%%%
bw2=zeros(ymax,xmax);
if(abs(theta-pi/2)<pi/4)
	x=[-xmax/2 xmax/2-1];
	y=rho/sin(theta)-x*cos(theta)/sin(theta);
	xindex=x+1+xmax/2;
	yindex=round(ymax/2-rho/sin(theta)+(xindex-1-xmax/2)*cos(theta)/sin(theta));
else
	y=[ymax/2-1 -ymax/2];
	x=rho/cos(theta)-y*sin(theta)/cos(theta);
	yindex=-y+ymax/2;
	xindex=round(rho/cos(theta)-(ymax/2-yindex)*sin(theta)/cos(theta)+xmax/2+1);
end
bw2=matline(bw2,yindex(1),xindex(1),yindex(2),xindex(2));
subplot(2,2,4);imshow(~bw2);
title(['逆hough変換結果 \rho=',num2str(rho),' \theta=',num2str(theta*180/pi)]);
subplot(2,2,2);imshow(~bw);title('現画像にオーバーラップ表示');
hold on
plot(xindex,yindex,'linewidth',5);
hold off
%drawnow;pause
%end
subplot(2,2,3);colormap hot;
