hough.m

%hough源函数
function [h,theta,rho]=hough(f,dtheta,drho)
%hough计算图像下f的Hough变换
%dtheta明确了在水平上的Hough变换的bins与theta坐标的空间。
%drho明确了在水平上Hough变换bins与rho坐标的空间。
%h是Hough变换矩阵。当nrho=2*ceil(norm(size(f))/drho)-1和ntheta=2*ceil(90/dtheta)时，它是nrho*ntheta。
%标记：如果90/dtheta不是一个整数，则实际角度是90/ceil（90/dtheta）。
%
%theta是一个包含h中每个列元素对应角度的ntheta元素向量。
%rho是一个包含h中每个元素对应的rho的值的nrho元素的向量。
%
%[h,theta,rho]=hough(f,dtheta,drho)使用dtheta=1和drho=1计算hough变换

if nargin<3
    drho=1;
end
if nargin<2
    dtheta=1;
end

f=double(f);
[M,N]=size(f);
theta=linspace(-90,0,ceil(90/dtheta)+1);
theta=[theta-fliplr(theta(2:end-1))];
ntheta=length(theta);
D=sqrt((M-1)^2+(N-1)^2);
q=ceil(D/drho);
nrho=2*q-1;
rho=linspace(-q*drho,q*drho,nrho);

[x,y,val]=find(f);
x=x-1;y-1;

%初始化输出
h=zeros(nrho,length(theta));

%防止过度记忆使用，每次执行1000个非零像素值
for k=1:ceil(length(val)/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,nthrta);
    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,1),1);
    
    theta_bin_index=repmat(1:ntheta,size(x_matrix,1),1);
    %实际的优势是：sparse函数的结构是一个sparse矩阵，当输入指标是重复时累加。这是我们想要Hough变换的一个动作。
    %我们想要输出一个满(非sparse)矩阵，然而，我们称函数full为sparse的输出
   h=h+full(sparse(rho_bin_index(:),theta_bin_index(:),...
       val_matrix(:),nrho,ntheta));
end