function lines = radonlines(varargin)
%RADONLINES Extract line segments based on Radon transform.
%
% 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.
% Author: HSW
% Date: 2015/4/21
% HARBIN INSTITUTE OF TECHNOLOGY
%
center = floor((0.5*size(varargin{1})));
centery = center(1);
centerx = center(2);
theta = varargin{2};
rho = varargin{3};
peaks = varargin{4};
type = varargin{5}; % 選擇畫直線還是畫線段
if isempty(peaks)disp('no peaks');return;
endif type == 1% 畫直線figure;imshow(varargin{1});hold on;for iter = 1:size(peaks,1)rho1 = rho(peaks(iter,1));theta1 = theta(peaks(iter,2));if rho1 <= 0 && theta1 <= 90% 直線在左下角theta2 = 90-theta1;slope = tan((theta2)*pi/180);y1 = centery - rho1*sin(theta1*pi/180);x1 = centerx + rho1*cos(theta1*pi/180);elseif rho1 < 0 && theta1 > 90% 直線在右下角theta2 = 270-theta1;theta1 = 180-theta1;slope = tan((theta2)*pi/180);y1 = centery - rho1*sin(theta1*pi/180);x1 = centerx - rho1*cos(theta1*pi/180);elseif rho1 > 0 && theta1 < 90% 直線在右上角theta2 = 90 - theta1;slope = tan((theta2)*pi/180);y1 = centery - rho1*sin(theta1*pi/180);x1 = centerx + rho1*cos(theta1*pi/180);else% 直線在左上角theta2 = 270-theta1;theta1 = 180-theta1;slope = tan((theta2)*pi/180);y1 = centery - rho1*sin(theta1*pi/180);x1 = centerx - rho1*cos(theta1*pi/180);endx =centerx-size(varargin{1},2):centerx+size(varargin{1},2);y = slope*(x-x1) + y1;plot(x,y,'g')end
elseif type == 2% 畫線段fillgap = varargin{6};minlength = varargin{7};delta = varargin{8};minlength_sq = minlength^2;fillgap_sq = fillgap^2;numlines = 0;[y,x] = find(varargin{1});nonzeropix = [x,y] - 1;lines = struct([]);for k = 1:size(peaks,1)[r,c] = radonpixels(nonzeropix,theta,rho,delta,peaks(k,:),center);if isempty(r)continue;end% Compute distance^2 between the point pairsxy = [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 thresholdfillgap_idx = find(dist_sq > fillgap_sq);idx = [0; fillgap_idx; size(xy,1)];for p = 1:length(idx) - 1p1 = xy(idx(p) + 1,:); % offset by 1 to convert to 1 based indexp2 = xy(idx(p + 1),:); % set the end (don't offset by one this time)linelength_sq = sum((p2-p1).^2);if linelength_sq >= minlength_sqnumlines = numlines + 1;lines(numlines).point1 = p1;lines(numlines).point2 = p2;lines(numlines).theta = theta(peaks(k,2));lines(numlines).rho = rho(peaks(k,1));endendend %for k = 1:size(peaks,1)
elseerror('type = 1 or type = 2');
end %if typeend %function radonlinesfunction [r,c] = radonpixels(nonzeropix,theta,rho,delta,peak,center)
x = nonzeropix(:,1);
y = nonzeropix(:,2);
centery = center(1);
centerx = center(2);
rho1 = rho(peak(1));
theta1 = theta(peak(2));
if rho1 <= 0 && theta1 <= 90% 直線在左下角theta2 = 90-theta1;slope = tan((theta2)*pi/180);y1 = centery - rho1*sin(theta1*pi/180);x1 = centerx + rho1*cos(theta1*pi/180);
elseif rho1 < 0 && theta1 > 90% 直線在右下角theta2 = 270-theta1;theta1 = 180-theta1;slope = tan((theta2)*pi/180);y1 = centery - rho1*sin(theta1*pi/180);x1 = centerx - rho1*cos(theta1*pi/180);
elseif rho1 > 0 && theta1 < 90% 直線在右上角theta2 = 90 - theta1;slope = tan((theta2)*pi/180);y1 = centery - rho1*sin(theta1*pi/180);x1 = centerx + rho1*cos(theta1*pi/180);
else% 直線在左上角theta2 = 270-theta1;theta1 = 180-theta1;slope = tan((theta2)*pi/180);y1 = centery - rho1*sin(theta1*pi/180);x1 = centerx - rho1*cos(theta1*pi/180);
end
idx = find(abs(slope*(x-x1) + y1 - y) <= delta); %進行直線擬合
r = y(idx) + 1;
c = x(idx) + 1;
[r,c] = reSortRadonPixels(r,c);
end% function radonpixelsfunction [r_new,c_new] = reSortRadonPixels(r,c)if isempty(r)r_new = r;c_new = c;return;
endr_range = max(r) - min(r);
c_range = max(c) - min(c);if r_range > c_rangesorting_order = [1,2];
elsesorting_order = [2,1];
end
[rc_new] = sortrows([r,c],sorting_order);
r_new = rc_new(:,1);
c_new = rc_new(:,2);
end % function reSortRadonPixels
function peaks = radonpeaks(varargin)
% RADONPEAKS Identify peaks in Radon 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].% odd 奇數
% '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.% Author: HSW
% HARBIN INSTITUTE OF TECHNOLOGY
[h, numpeaks, threshold, nhood] = parseInputs(varargin{:});
% h: radon 變換的輸出
% numpeaks: 峰值的個數
% threshold: 峰值的最小值, 默認為0.5*max(H(:))
% nhood: 包含兩個奇數的數組[M,N], 當峰值識別出來后,設置為0 % initialize the loop variables
done = false;
hnew = h;
nhood_center = (nhood-1)/2;% 抑制塊的中心位置,例如nhood = [5,5], 則nhood_center = [2,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) >= thresholdpeaks = [peaks; [p q]]; %#ok<AGROW> % 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; %設置為0 done = size(peaks,1) == numpeaks;elsedone = true;end
end%-----------------------------------------------------------------------------
function [h, numpeaks, threshold, nhood] = parseInputs(varargin)narginchk(1,6); % 參數個數小于1 大于6報錯h = varargin{1};
% validateattributes(h, {'double'}, {'real', '2d', 'nonsparse', 'nonempty',...
% 'finite', 'integer'}, ...
% mfilename, 'H', 1);
% hough變換的中h中的取值必然為非負整數,但是,radon變換中只能保證為非負的validateattributes(h, {'double'}, {'real', '2d', 'nonsparse', 'nonempty',...'finite'}, ...mfilename, 'H', 1);numpeaks = 1; % set default value for numpeaks峰值的默認值為1idx = 2;
if nargin > 1if ~ischar(varargin{2})numpeaks = varargin{2};validateattributes(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 pairsdone = false;while ~doneinput = varargin{idx};inputStr = validatestring(input, validStrings,mfilename,'PARAM',idx);idx = idx+1; %advance index to point to the VAL portion of the input if idx > narginerror(message('images:houghpeaks:valForhoughpeaksMissing', inputStr))endswitch inputStrcase 'Threshold'threshold = varargin{idx};validateattributes(threshold, {'double'}, {'real', 'scalar','nonnan' ...'nonnegative'}, mfilename, inputStr, idx);case 'NHoodSize'nhood = varargin{idx};validateattributes(nhood, {'double'}, {'real', 'vector', ...'finite','integer','positive','odd'},...mfilename, inputStr, idx);if (any(size(nhood) ~= [1,2]))error(message('images:radonpeaks:invalidNHoodSize', inputStr))endif (any(nhood > size(h)))error(message('images:radonpeaks:tooBigNHoodSize', inputStr))end otherwise%should never get hereerror(message('images:radonpeaks:internalError'))endif idx >= nargindone = true;endidx=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.確保nhood為奇數
endif isempty(threshold)threshold = 0.5 * max(h(:)); %設置默認值
end
