%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