function peaks = houghpeaks(varargin)
%HOUGHPEAKS Identify peaks in the Hough transform.
% PEAKS = HOUGHPEAKS(H,NUMPEAKS) locates peaks in the Hough
% transform matrix, H, generated by the HOUGH function. NUMPEAKS
% specifies the maximum number of peaks to identify. PEAKS is
% a Q-by-2 matrix, where Q can range from 0 to NUMPEAKS. Q holds
% the row and column coordinates of the peaks. If NUMPEAKS is
% omitted, it defaults to 1.
%
% PEAKS = HOUGHPEAKS(...,PARAM1,VAL1,PARAM2,VAL2) sets various
% parameters. Parameter names can be abbreviated, and case
% does not matter. Each string parameter is followed by a value
% as indicated below:
%
% 'Threshold' Nonnegative scalar.
% Values of H below 'Threshold' will not be considered
% to be peaks. Threshold can vary from 0 to Inf.
%
% Default: 0.5*max(H(:))
%
% 'NHoodSize' Two-element vector of positive odd integers: [M N].
% 'NHoodSize' specifies the size of the suppression
% neighborhood. This is the neighborhood around each
% peak that is set to zero after the peak is identified.
%
% Default: smallest odd values greater than or equal to
% size(H)/50.
%
% Class Support
% -------------
% H is the output of the HOUGH function. NUMPEAKS is a positive
% integer scalar.
%
% Example
% -------
% Locate and display two peaks in the Hough transform of the
% rotated circuit.tif image.
%
% I = imread('circuit.tif');
% BW = edge(imrotate(I,50,'crop'),'canny');
% [H,T,R] = hough(BW);
% P = houghpeaks(H,2);
% imshow(H,[],'XData',T,'YData',R,'InitialMagnification','fit');
% xlabel('\theta'), ylabel('\rho');
% axis on, axis normal, hold on;
% plot(T(P(:,2)),R(P(:,1)),'s','color','white');
%
% See also HOUGH and HOUGHLINES.
% Copyright 1993-2003 The MathWorks, Inc.
% $Revision: 1.1.8.2 $ $Date: 2004/10/20 17:54:16 $
% References:
% Rafael C. Gonzalez, Richard E. Woods, Steven L. Eddins, "Digital
% Image Processing Using MATLAB", Prentice Hall, 2004
[h, numpeaks, threshold, nhood] = parseInputs(varargin{:});
% initialize the loop variables
done = false;
hnew = h;
nhood_center = (nhood-1)/2;
peaks = [];
while ~done
[dummy max_idx] = max(hnew(:)); %#ok
[p, q] = ind2sub(size(hnew), max_idx);
p = p(1); q = q(1);
if hnew(p, q) >= threshold
peaks = [peaks; [p q]]; % add the peak to the list
% Suppress this maximum and its close neighbors.
p1 = p - nhood_center(1); p2 = p + nhood_center(1);
q1 = q - nhood_center(2); q2 = q + nhood_center(2);
% Throw away neighbor coordinates that are out of bounds in
% the rho direction.
[qq, pp] = meshgrid(q1:q2, max(p1,1):min(p2,size(h,1)));
pp = pp(:); qq = qq(:);
% For coordinates that are out of bounds in the theta
% direction, we want to consider that H is antisymmetric
% along the rho axis for theta = +/- 90 degrees.
theta_too_low = find(qq < 1);
qq(theta_too_low) = size(h, 2) + qq(theta_too_low);
pp(theta_too_low) = size(h, 1) - pp(theta_too_low) + 1;
theta_too_high = find(qq > size(h, 2));
qq(theta_too_high) = qq(theta_too_high) - size(h, 2);
pp(theta_too_high) = size(h, 1) - pp(theta_too_high) + 1;
% Convert to linear indices to zero out all the values.
hnew(sub2ind(size(hnew), pp, qq)) = 0;
done = size(peaks,1) == numpeaks;
else
done = true;
end
end
%-----------------------------------------------------------------------------
function [h, numpeaks, threshold, nhood] = parseInputs(varargin)
iptchecknargin(1,6,nargin,mfilename);
h = varargin{1};
iptcheckinput(h, {'double'}, {'real', '2d', 'nonsparse', 'nonempty',...
'finite', 'integer'}, ...
mfilename, 'H', 1);
numpeaks = 1; % set default value for numpeaks
idx = 2;
if nargin > 1
if ~ischar(varargin{2})
numpeaks = varargin{2};
iptcheckinput(numpeaks, {'double'}, {'real', 'scalar', 'integer', ...
'positive', 'nonempty'}, mfilename, 'NUMPEAKS', 2);
idx = 3;
end
end
% Initialize to empty
nhood = [];
threshold = [];
% Process parameter-value pairs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
validStrings = {'Threshold','NHoodSize'};
if nargin > idx-1 % we have parameter/value pairs
done = false;
while ~done
input = varargin{idx};
inputStr = iptcheckstrs(input, validStrings,mfilename,'PARAM',idx);
idx = idx+1; %advance index to point to the VAL portion of the input
if idx > nargin
eid = sprintf('Images:%s:valFor%sMissing', mfilename, inputStr);
msg = sprintf('Parameter ''%s'' must be followed by a value.', inputStr);
error(eid,'%s', msg);
end
switch inputStr
case 'Threshold'
threshold = varargin{idx};
iptcheckinput(threshold, {'double'}, {'real', 'scalar','nonnan' ...
'nonnegative'}, mfilename, inputStr, idx);
case 'NHoodSize'
nhood = varargin{idx};
iptcheckinput(nhood, {'double'}, {'real', 'vector', ...
'finite','integer','positive','odd'},...
mfilename, inputStr, idx);
if (any(size(nhood) ~= [1,2]))
eid = sprintf('Images:%s:invalidNHoodSize', mfilename);
msg = sprintf(['Value of parameter,''%s'', must be a two '...
'element row vector.'], inputStr);
error(eid,'%s', msg);
end
if (any(nhood > size(h)))
eid = sprintf('Images:%s:tooBigNHoodSize', mfilename);
msg = sprintf(['Dimensions specified by ''%s'', '...
'must be smaller than '...
'size of Hough matrix H.'], inputStr);
error(eid,'%s', msg);
end
otherwise
eid = sprintf('Images:%s:internalError', mfilename);
msg = 'Unknown input string.'; %should never get here
error(eid,'%s', msg);
end
if idx >= nargin
done = true;
end
idx=idx+1;
end
end
% Set the defaults if necessary
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if isempty(nhood)
nhood = size(h)/50;
nhood = max(2*ceil(nhood/2) + 1, 1); % Make sure the nhood size is odd.
end
if isempty(threshold)
threshold = 0.5 * max(h(:));
end
function lines = houghlines(varargin)
%HOUGHLINES Extract line segments based on Hough transform.
% LINES = HOUGHLINES(BW, THETA, RHO, PEAKS) extracts line segments
% in the image BW associated with particular bins in a Hough
% transform. THETA and RHO are vectors returned by function HOUGH.
% Matrix PEAKS, which is returned by function HOUGHPEAKS,
% contains the row and column coordinates of the Hough transform
% bins to use in searching for line segments. HOUGHLINES returns
% LINES structure array whose length equals the number of merged
% line segments found. Each element of the structure array has
% these fields:
%
% point1 End-point of the line segment; two-element vector
% point2 End-point of the line segment; two-element vector
% theta Angle (in degrees) of the Hough transform bin
% rho Rho-axis position of the Hough transform bin
%
% The end-point vectors contain [X, Y] coordinates.
%
% LINES = HOUGHLINES(...,PARAM1,VAL1,PARAM2,VAL2) sets various
% parameters. Parameter names can be abbreviated, and case
% does not matter. Each string parameter is followed by a value
% as indicated below:
%
% 'FillGap' Positive real scalar.
% When HOUGHLINES finds two line segments associated
% with the same Hough transform bin that are separated
% by less than 'FillGap' distance, HOUGHLINES merges
% them into a single line segment.
%
% Default: 20
%
% 'MinLength' Positive real scalar.
% Merged line segments shorter than 'MinLength'
% are discarded.
%
% Default: 40
%
% Class Support
% -------------
% BW can be logical or numeric and it must be real, 2-D, and nonsparse.
%
% Example
% -------
% Search for line segments corresponding to five peaks in the Hough
% transform of the rotated circuit.tif image. Additionally, highlight
% the longest segment.
%
% I = imread('circuit.tif');
% rotI = imrotate(I,33,'crop');
% BW = edge(rotI,'canny');
% [H,T,R] = hough(BW);
% imshow(H,[],'XData',T,'YData',R,'InitialMagnification','fit');
% xlabel('\theta'), ylabel('\rho');
% axis on, axis normal, hold on;
% P = houghpeaks(H,5,'threshold',ceil(0.3*max(H(:))));
% x = T(P(:,2)); y = R(P(:,1));
% plot(x,y,'s','color','white');
%
% % Find lines and plot them
% lines = houghlines(BW,T,R,P,'FillGap',5,'MinLength',7);
% figure, imshow(rotI), hold on
% max_len = 0;
% for k = 1:length(lines)
% xy = [lines(k).point1; lines(k).point2];
% plot(xy(:,1),xy(:,2),'LineWidth',2,'Color','green');
%
% % plot beginnings and ends of lines
% plot(xy(1,1),xy(1,2),'x','LineWidth',2,'Color','yellow');
% plot(xy(2,1),xy(2,2),'x','LineWidth',2,'Color','red');
%
% % determine the endpoints of the longest line segment
% len = norm(lines(k).point1 - lines(k).point2);
% if ( len > max_len)
% max_len = len;
% xy_long = xy;
% end
% end
%
% % highlight the longest line segment
% plot(xy_long(:,1),xy_long(:,2),'LineWidth',2,'Color','cyan');
%
% See also HOUGH and HOUGHPEAKS.
% Copyright 1993-2003 The MathWorks, Inc.
% $Revision: 1.1.8.2 $ $Date: 2005/12/12 23:20:08 $
% References:
% Rafael C. Gonzalez, Richard E. Woods, Steven L. Eddins, "Digital
% Image Processing Using MATLAB", Prentice Hall, 2003
[nonzeropix,theta,rho,peaks,fillgap,minlength] = parseInputs(varargin{:});
minlength_sq = minlength^2;
fillgap_sq = fillgap^2;
numlines = 0; lines = struct;
for k = 1:size(peaks,1)
% Get all pixels associated with Hough transform cell.
[r, c] = houghpixels(nonzeropix, theta, rho, peaks(k,:));
if isempty(r)
continue
end
% Compute distance^2 between the point pairs
xy = [c r]; % x,y pairs in coordinate system with the origin at (1,1)
diff_xy_sq = diff(xy,1,1).^2;
dist_sq = sum(diff_xy_sq,2);
% Find the gaps larger than the threshold.
fillgap_idx = find(dist_sq > fillgap_sq);
idx = [0; fillgap_idx; size(xy,1)];
for p = 1:length(idx) - 1
p1 = xy(idx(p) + 1,:); % offset by 1 to convert to 1 based index
p2 = xy(idx(p + 1),:); % set the end (don't offset by one this time)
linelength_sq = sum((p2-p1).^2);
if linelength_sq >= minlength_sq
numlines = numlines + 1;
lines(numlines).point1 = p1;
lines(numlines).point2 = p2;
lines(numlines).theta = theta(peaks(k,2));
lines(numlines).rho = rho(peaks(k,1));
end
end
end
%-----------------------------------------------------------------------------
function [r, c] = houghpixels(nonzeropix, theta, rho, peak)
%HOUGHPIXELS Compute image pixels belonging to Hough transform bin.
% [R, C] = HOUGHPIXELS(NONZEROPIX, THETA, RHO, PEAK) computes the
% row-column indices (R, C) for nonzero pixels NONZEROPIX that map
% to a particular Hough transform bin, PEAK which is a two element
% vector [RBIN CBIN]. RBIN and CBIN are scalars indicating the
% row-column bin location in the Hough transform matrix returned by
% function HOUGH. THETA and RHO are the second and third output
% arguments from the HOUGH function.
x = nonzeropix(:,1);
y = nonzeropix(:,2);
theta_c = theta(peak(2)) * pi / 180;
rho_xy = x*cos(theta_c) + y*sin(theta_c);
nrho = length(rho);
slope = (nrho - 1)/(rho(end) - rho(1));
rho_bin_index = round(slope*(rho_xy - rho(1)) + 1);
idx = find(rho_bin_index == peak(1));
r = y(idx) + 1; c = x(idx) + 1;
%-----------------------------------------------------------------------------
function [nonzeropix,theta,rho,peaks,fillgap,minlength] = ...
parseInputs(varargin)
iptchecknargin(1,8,nargin,mfilename);
idx = 1;
bw = varargin{idx};
iptcheckinput(bw, {'numeric','logical'},...
{'real', '2d', 'nonsparse', 'nonempty'}, ...
mfilename, 'BW', idx);
idx = idx+1;
theta = varargin{idx};
iptcheckinput(theta, {'double'}, {'real','vector','finite',...
'nonsparse','nonempty'}, ...
mfilename, 'THETA', idx);
idx = idx+1;
rho = varargin{idx};
iptcheckinput(rho, {'double'}, {'real','vector','finite',...
'nonsparse','nonempty'}, ...
mfilename, 'RHO', idx);
idx = idx+1;
peaks = varargin{idx};
iptcheckinput(peaks, {'double'}, {'real','2d','nonsparse','integer'}, ...
mfilename, 'PEAKS', idx);
if size(peaks,2) ~= 2
eid = sprintf('Images:%s:invalidPEAKS', mfilename);
msg = sprintf('PEAKS must be a Q-by-2 matrix');
error(eid,'%s',msg);
end
% Set the defaults
fillgap = 20;
minlength = 40;
% Process parameter-value pairs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
validStrings = {'FillGap','MinLength'};
idx = idx+1;
if nargin > idx-1 % we have parameter/value pairs
done = false;
while ~done
input = varargin{idx};
inputStr = iptcheckstrs(input, validStrings,mfilename,'PARAM',idx);
idx = idx+1; %advance index to point to the VAL portion of the input
if idx > nargin
eid = sprintf('Images:%s:valFor%sMissing', mfilename, inputStr);
msg = sprintf('Parameter ''%s'' must be followed by a value.', inputStr);
error(eid,'%s',msg);
end
switch inputStr
case 'FillGap'
fillgap = varargin{idx};
iptcheckinput(fillgap, {'double'}, {'finite','real', 'scalar', ...
'positive'}, mfilename, inputStr, idx);
case 'MinLength'
minlength = varargin{idx};
iptcheckinput(minlength, {'double'}, {'finite','real', 'scalar', ...
'positive'}, mfilename, inputStr, idx);
otherwise
eid = sprintf('Images:%s:internalError', mfilename);
msg = 'Unknown input string.'; %should never get here
error(eid,'%s',msg);
end
if idx >= nargin
done = true;
end
idx=idx+1;
end
end
% Compute the required parameters
[y, x] = find(bw);
nonzeropix = [x, y] - 1;
简化整合为一个hough变换函数