houghpeaks

%houghpeaks源函数
function [r,c,hnew]=houghpeaks(h,numpeaks,threshold,nhood)
%houghpeaks在Hough中检测峰值.
%[r,c,hnew]=houghpeaks(h,numpeaks,threshold,nhood)在Hough矩阵h中检测峰值.
%numpeaks明确地寻找最大值的位置。
%在threshold下h的值不会被咯频率为峰值。
%nhood是一个两个元素的向量明确了抑制领域的大小。
%当峰值被确立了，每个峰值的领域设为0。
%nhood的元素必须是正的，奇数的整数。
%r和c是确立峰值的行坐标和列坐标。
%
%若nhood被忽略，则默认为大于size(h)/50的最小奇数。
%若threshold被忽略，则默认为0.5*max(h(:))。
%若numpeaks被忽略，则默认为1.

if nargin<4
    nhood=size(h)/50;
    %确保领域是奇的。
    nhood=max(2*ceil(nhood/2)+1,1);
end
if nargin<3
    threshold=0.5*max(h(:));
end
if nargin<2
    numpeaks=1;
end

done=false;
hnew=h;r=[];c=[];
while~done
    [p,q]=find(hnew==max(hnew(:)));
    p=p(1);q=q(1);
    if hnew(p,q)>=threshold
        r(end+1)=p;c(end+1)=q;
        
        %印制最大值和他的领域。
        p1=p-(nhood(1)-1)/2;p2=p+(nhood(1)-1)/2;
        q1=q-(nhood(2)-1)/2;q2=q+(nhood(2)-1)/2;
        [pp,qq]=ndgrid(p1:p2,q1:q2);
        
        %当在rho方向超出边界，扔掉领域坐标。
        badrho=find((pp<1)|(pp>size(h,1)));
        pp(badrho)=[];qq(badrho)=[];
        %为了在theta方向超出边界的坐标，我们要考虑h在rho坐标是反对称的，因为theta=+/-90.        
        theta_too_low=find(qq<1);
        qq(theta_too_low)=size(h,2)+qq(theta_too_low);
        pp(theta_too_low)=size(h,1)-qq(theta_too_low)+1;
        theta_too_low=find(qq>size(h,2));
        qq(theta_too_low)=qq(theta_too_high)-size(h,2);
        pp(theta_too_low)=size(h,1)-qq(theta_too_high)+1;
        
        %把所有值转线性指标为0.
        hnew(sub2ind(size(hnew),pp,qq))=0;
        done=length(r)==numpeaks;
    else
        done=true;
    end
end