clc;
clear;

RGB = imread('test.jpg');%jaynes-thesis
I=rgb2gray(RGB);  % 圖片用的是灰度圖像，
[x,y,z]=size(I); 
BW=edge(I,'canny');

figure;
imshow(I);
title('原圖');

figure;
imshow(BW);
title('邊緣檢測(cè)圖像');

rho_max=floor(sqrt(x^2+y^2))+1;  %由原圖數(shù)組坐標(biāo)算出ρ最大值，并取整數(shù)部分加1
                                 %此值作為ρ，θ坐標(biāo)系ρ最大值
accarray=zeros(rho_max,180);     %定義ρ，θ坐標(biāo)系的數(shù)組，初值為0。
%θ的最大值，180度

Theta=[0:pi/180:pi];             %定義θ數(shù)組，確定θ取值范圍

for n=1:x,
   for m=1:y
      if BW(n,m)==1
         for k=1:180
            %將θ值代入hough變換方程，求ρ值
            rho=(m*cos(Theta(k)))+(n*sin(Theta(k)));
            %將ρ值與ρ最大值的和的一半作為ρ的坐標(biāo)值（數(shù)組坐標(biāo)），這樣做是為了防止ρ值出現(xiàn)負(fù)數(shù)
            rho_int=round(rho/2+rho_max/2);
            %在ρθ坐標(biāo)（數(shù)組）中標(biāo)識(shí)點(diǎn)，即計(jì)數(shù)累加
            accarray(rho_int,k)=accarray(rho_int,k)+1;
         end
      end
   end
end

%figure;colormap gray;
%imagesc(accarray);title('hough變換后的圖')
%xlabel('\theta'), ylabel('\rho');

%=====下面程序的顯示效果沒(méi)上面好=====%
%accarray=uint8(accarray);                %轉(zhuǎn)換后會(huì)丟數(shù)據(jù)
%figure;imshow(accarray);title('hough變換后的圖')
%xlabel('\theta'), ylabel('\rho');
%axis on, axis normal, hold on;

%=======利用hough變換提取直線======%
%尋找0.25%的圖像大小個(gè)數(shù)像素以上的直線在hough變換后形成的點(diǎn)
K=1;                    %存儲(chǔ)數(shù)組計(jì)數(shù)器
for rho_n=1:rho_max     %在hough變換后的數(shù)組中搜索
    for theta_m=1:180
        if accarray(rho_n,theta_m)>=20 & rho_n<rho_max&accarray(rho_n,theta_m)>=accarray(rho_n+1,theta_m) & rho_n>1&accarray(rho_n,theta_m)>=accarray(rho_n-1,theta_m) & theta_m>1&accarray(rho_n,theta_m)>=accarray(rho_n,theta_m-1) & theta_m<179&accarray(rho_n,theta_m)>=accarray(rho_n,theta_m+1)%設(shè)定直線的最小值。
            case_accarray_n(K)=rho_n;     %存儲(chǔ)搜索出的數(shù)組下標(biāo)
            case_accarray_m(K)=theta_m;
            K=K+1;
        end
    end
end
%把這些點(diǎn)構(gòu)成的直線提取出來(lái),輸出圖像數(shù)組為I_out
I_out=zeros(x,y);
I_jiao_class=zeros(x,y);
for n=1:x,
   for m=1:y
      if BW(n,m)==1
         for k=1:180
            rho=(m*cos(Theta(k)))+(n*sin(Theta(k)));
            rho_int=round(rho/2+rho_max/2);
            %如果正在計(jì)算的點(diǎn)屬于0.25%圖像大小像素以上點(diǎn)，則把它提取出來(lái)
            for a=1:K-1
                if rho_int==case_accarray_n(a)&k==case_accarray_m(a)%%%==gai==%%% k==case_accarray_m(a)&rho_int==case_accarray_n(a)
                    I_out(n,m)=BW(n,m);  
                    I_jiao_class(n,m)=k;
                end
            end
         end
      end
   end
end

figure;
imshow(I_out);
title('利用經(jīng)典hough變換提取的圖像');

%===============New    Code==================%
%disp(case_accarray_n);
%disp(case_accarray_m);
%dignant_k=zeros(1,K);
%dignant_b=zeros(1,K);
for k=1:K-1
    if case_accarray_m~=0
        dignant_k(k)=-atan(case_accarray_m(k));
        dignant_b(k)=case_accarray_n(k)/sin(case_accarray_m(k));
    else
        dignant_k(k)=Inf;
        dignant_b(k)=case_accarray_n(k);
    end
end
disp(dignant_k);
disp(dignant_b);