function [h,theta,rho]=hough_trans(f,th0,p0)

f=double(f);

[M,N]=size(f);
theta=linspace(-90,0,ceil(90/th0)+1);%產生theta坐標系
theta=[theta -fliplr(theta(2:end-1))];
ntheta=length(theta);

D=sqrt((M-1)^2+(N-1)^2);
q=ceil(D/p0);
nrho=2*q-1;
rho=linspace(-q*p0,q*p0,nrho);%建立rho坐標
[x,y,val]=find(f);%找出邊緣的坐標，及其相應灰度值1
x=x-1;y=y-1;%坐標平移


h=zeros(nrho,ntheta);%初始化hough圖矩陣

for k=1:ceil(length(x)/1000)%每1000個處理一次
    first=(k-1)*1000+1;
    last=min(first+999,length(x));
    x_matrix=repmat(x(first:last),1,ntheta);
    y_matrix=repmat(y(first:last),1,ntheta);
    val_matrix=repmat(val(first:last),1,ntheta);%不必要
     
    theta_matrix=repmat(theta,size(x_matrix,1),1)*pi/180;
    
    rho_matrix=x_matrix.*cos(theta_matrix)+y_matrix.*sin(theta_matrix);
    slope=(nrho-1)/(rho(end)-rho(1));
    rho_bin_index=round(slope*(rho_matrix-rho(1))+1);
    theta_bin_index=repmat(1:ntheta,size(x_matrix,1),1);
    h=h+full(sparse(rho_bin_index(:),theta_bin_index(:),val_matrix(:),nrho,ntheta));
end
% for i=1:length(x)循環太慢
% for j=1:length(y)
%      
% for k=1:length(theta)
% p=round(x(i)*cos(theta(k)*pi/180)+y(j)*sin(theta(k)*pi/180));
% m=round((theta(k)+90)/th0)+1;n=round((p+q)/p0)+1;
% h(n,m)=h(n,m)+1;
%  
% end;
%     end;
% end;