www.gusucode.com > 基于matlab程序hough变换检测圆源码 > 基于matlab程序hough变换检测圆源码/hough检测圆/CircularHough_Grd.m
function [accum, varargout] = CircularHough_Grd(img, radrange, varargin) %Detect circular shapes in a grayscale image. Resolve their center %positions and radii. % % [accum, circen, cirrad, dbg_LMmask] = CircularHough_Grd( % img, radrange, grdthres, fltr4LM_R, multirad, fltr4accum) % Circular Hough transform based on the gradient field of an image. % NOTE: Operates on grayscale images, NOT B/W bitmaps. % NO loops in the implementation of Circular Hough transform, % which means faster operation but at the same time larger % memory consumption. % %%%%%%%% INPUT: (img, radrange, grdthres, fltr4LM_R, multirad, fltr4accum) % % img: A 2-D grayscale image (NO B/W bitmap) % % radrange: The possible minimum and maximum radii of the circles % to be searched, in the format of % [minimum_radius , maximum_radius] (unit: pixels) % **NOTE**: A smaller range saves computational time and % memory. % % grdthres: (Optional, default is 10, must be non-negative) % The algorithm is based on the gradient field of the % input image. A thresholding on the gradient magnitude % is performed before the voting process of the Circular % Hough transform to remove the 'uniform intensity' % (sort-of) image background from the voting process. % In other words, pixels with gradient magnitudes smaller % than 'grdthres' are NOT considered in the computation. % **NOTE**: The default parameter value is chosen for % images with a maximum intensity close to 255. For cases % with dramatically different maximum intensities, e.g. % 10-bit bitmaps in stead of the assumed 8-bit, the default % value can NOT be used. A value of 4% to 10% of the maximum % intensity may work for general cases. % % fltr4LM_R: (Optional, default is 8, minimum is 3) % The radius of the filter used in the search of local % maxima in the accumulation array. To detect circles whose % shapes are less perfect, the radius of the filter needs % to be set larger. % % multirad: (Optional, default is 0.5) % In case of concentric circles, multiple radii may be % detected corresponding to a single center position. This % argument sets the tolerance of picking up the likely % radii values. It ranges from 0.1 to 1, where 0.1 % corresponds to the largest tolerance, meaning more radii % values will be detected, and 1 corresponds to the smallest % tolerance, in which case only the "principal" radius will % be picked up. % % fltr4accum: (Optional. A default filter will be used if not given) % Filter used to smooth the accumulation array. Depending % on the image and the parameter settings, the accumulation % array built has different noise level and noise pattern % (e.g. noise frequencies). The filter should be set to an % appropriately size such that it's able to suppress the % dominant noise frequency. % %%%%%%%% OUTPUT: [accum, circen, cirrad, dbg_LMmask] % % accum: The result accumulation array from the Circular Hough % transform. The accumulation array has the same dimension % as the input image. % % circen: (Optional) % Center positions of the circles detected. Is a N-by-2 % matrix with each row contains the (x, y) positions % of a circle. For concentric circles (with the same center % position), say k of them, the same center position will % appear k times in the matrix. % % cirrad: (Optional) % Estimated radii of the circles detected. Is a N-by-1 % column vector with a one-to-one correspondance to the % output 'circen'. A value 0 for the radius indicates a % failed detection of the circle's radius. % % dbg_LMmask: (Optional, for debugging purpose) % Mask from the search of local maxima in the accumulation % array. % %%%%%%%%% EXAMPLE #0: % rawimg = imread('TestImg_CHT_a2.bmp'); % tic; % [accum, circen, cirrad] = CircularHough_Grd(rawimg, [15 60]); % toc; % figure(1); imagesc(accum); axis image; % title('Accumulation Array from Circular Hough Transform'); % figure(2); imagesc(rawimg); colormap('gray'); axis image; % hold on; % plot(circen(:,1), circen(:,2), 'r+'); % for k = 1 : size(circen, 1), % DrawCircle(circen(k,1), circen(k,2), cirrad(k), 32, 'b-'); % end % hold off; % title(['Raw Image with Circles Detected ', ... % '(center positions and radii marked)']); % figure(3); surf(accum, 'EdgeColor', 'none'); axis ij; % title('3-D View of the Accumulation Array'); % % COMMENTS ON EXAMPLE #0: % Kind of an easy case to handle. To detect circles in the image whose % radii range from 15 to 60. Default values for arguments 'grdthres', % 'fltr4LM_R', 'multirad' and 'fltr4accum' are used. % %%%%%%%%% EXAMPLE #1: % rawimg = imread('TestImg_CHT_a3.bmp'); % tic; % [accum, circen, cirrad] = CircularHough_Grd(rawimg, [15 60], 10, 20); % toc; % figure(1); imagesc(accum); axis image; % title('Accumulation Array from Circular Hough Transform'); % figure(2); imagesc(rawimg); colormap('gray'); axis image; % hold on; % plot(circen(:,1), circen(:,2), 'r+'); % for k = 1 : size(circen, 1), % DrawCircle(circen(k,1), circen(k,2), cirrad(k), 32, 'b-'); % end % hold off; % title(['Raw Image with Circles Detected ', ... % '(center positions and radii marked)']); % figure(3); surf(accum, 'EdgeColor', 'none'); axis ij; % title('3-D View of the Accumulation Array'); % % COMMENTS ON EXAMPLE #1: % The shapes in the raw image are not very good circles. As a result, % the profile of the peaks in the accumulation array are kind of % 'stumpy', which can be seen clearly from the 3-D view of the % accumulation array. (As a comparison, please see the sharp peaks in % the accumulation array in example #0) To extract the peak positions % nicely, a value of 20 (default is 8) is used for argument 'fltr4LM_R', % which is the radius of the filter used in the search of peaks. % %%%%%%%%% EXAMPLE #2: % rawimg = imread('TestImg_CHT_b3.bmp'); % fltr4img = [1 1 1 1 1; 1 2 2 2 1; 1 2 4 2 1; 1 2 2 2 1; 1 1 1 1 1]; % fltr4img = fltr4img / sum(fltr4img(:)); % imgfltrd = filter2( fltr4img , rawimg ); % tic; % [accum, circen, cirrad] = CircularHough_Grd(imgfltrd, [15 80], 8, 10); % toc; % figure(1); imagesc(accum); axis image; % title('Accumulation Array from Circular Hough Transform'); % figure(2); imagesc(rawimg); colormap('gray'); axis image; % hold on; % plot(circen(:,1), circen(:,2), 'r+'); % for k = 1 : size(circen, 1), % DrawCircle(circen(k,1), circen(k,2), cirrad(k), 32, 'b-'); % end % hold off; % title(['Raw Image with Circles Detected ', ... % '(center positions and radii marked)']); % % COMMENTS ON EXAMPLE #2: % The circles in the raw image have small scale irregularities along % the edges, which could lead to an accumulation array that is bad for % local maxima detection. A 5-by-5 filter is used to smooth out the % small scale irregularities. A blurred image is actually good for the % algorithm implemented here which is based on the image's gradient % field. % %%%%%%%%% EXAMPLE #3: % rawimg = imread('TestImg_CHT_c3.bmp'); % fltr4img = [1 1 1 1 1; 1 2 2 2 1; 1 2 4 2 1; 1 2 2 2 1; 1 1 1 1 1]; % fltr4img = fltr4img / sum(fltr4img(:)); % imgfltrd = filter2( fltr4img , rawimg ); % tic; % [accum, circen, cirrad] = ... % CircularHough_Grd(imgfltrd, [15 105], 8, 10, 0.7); % toc; % figure(1); imagesc(accum); axis image; % figure(2); imagesc(rawimg); colormap('gray'); axis image; % hold on; % plot(circen(:,1), circen(:,2), 'r+'); % for k = 1 : size(circen, 1), % DrawCircle(circen(k,1), circen(k,2), cirrad(k), 32, 'b-'); % end % hold off; % title(['Raw Image with Circles Detected ', ... % '(center positions and radii marked)']); % % COMMENTS ON EXAMPLE #3: % Similar to example #2, a filtering before circle detection works for % noisy image too. 'multirad' is set to 0.7 to eliminate the false % detections of the circles' radii. % %%%%%%%%% BUG REPORT: % This is a beta version. Please send your bug reports, comments and % suggestions to pengtao@glue.umd.edu . Thanks. % % %%%%%%%%% INTERNAL PARAMETERS: % The INPUT arguments are just part of the parameters that are used by % the circle detection algorithm implemented here. Variables in the code % with a prefix 'prm_' in the name are the parameters that control the % judging criteria and the behavior of the algorithm. Default values for % these parameters can hardly work for all circumstances. Therefore, at % occasions, the values of these INTERNAL PARAMETERS (parameters that % are NOT exposed as input arguments) need to be fine-tuned to make % the circle detection work as expected. % The following example shows how changing an internal parameter could % influence the detection result. % 1. Change the value of the internal parameter 'prm_LM_LoBndRa' to 0.4 % (default is 0.2) % 2. Run the following matlab code: % fltr4accum = [1 2 1; 2 6 2; 1 2 1]; % fltr4accum = fltr4accum / sum(fltr4accum(:)); % rawimg = imread('Frame_0_0022_portion.jpg'); % tic; % [accum, circen] = CircularHough_Grd(rawimg, ... % [4 14], 10, 4, 0.5, fltr4accum); % toc; % figure(1); imagesc(accum); axis image; % title('Accumulation Array from Circular Hough Transform'); % figure(2); imagesc(rawimg); colormap('gray'); axis image; % hold on; plot(circen(:,1), circen(:,2), 'r+'); hold off; % title('Raw Image with Circles Detected (center positions marked)'); % 3. See how different values of the parameter 'prm_LM_LoBndRa' could % influence the result. % Author: Tao Peng % Department of Mechanical Engineering % University of Maryland, College Park, Maryland 20742, USA % pengtao@glue.umd.edu % Version: Beta Revision: Mar. 07, 2007 %%%%%%%% Arguments and parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Validation of arguments if ndims(img) ~= 2 || ~isnumeric(img), error('CircularHough_Grd: ''img'' has to be 2 dimensional'); end if ~all(size(img) >= 32), error('CircularHough_Grd: ''img'' has to be larger than 32-by-32'); end if numel(radrange) ~= 2 || ~isnumeric(radrange), error(['CircularHough_Grd: ''radrange'' has to be ', ... 'a two-element vector']); end prm_r_range = sort(max( [0,0;radrange(1),radrange(2)] )); % Parameters (default values) prm_grdthres = 10; prm_fltrLM_R = 8; prm_multirad = 0.5; func_compu_cen = true; func_compu_radii = true; % Validation of arguments vap_grdthres = 1; if nargin > (1 + vap_grdthres), if isnumeric(varargin{vap_grdthres}) && ... varargin{vap_grdthres}(1) >= 0, prm_grdthres = varargin{vap_grdthres}(1); else error(['CircularHough_Grd: ''grdthres'' has to be ', ... 'a non-negative number']); end end vap_fltr4LM = 2; % filter for the search of local maxima if nargin > (1 + vap_fltr4LM), if isnumeric(varargin{vap_fltr4LM}) && varargin{vap_fltr4LM}(1) >= 3, prm_fltrLM_R = varargin{vap_fltr4LM}(1); else error(['CircularHough_Grd: ''fltr4LM_R'' has to be ', ... 'larger than or equal to 3']); end end vap_multirad = 3; if nargin > (1 + vap_multirad), if isnumeric(varargin{vap_multirad}) && ... varargin{vap_multirad}(1) >= 0.1 && ... varargin{vap_multirad}(1) <= 1, prm_multirad = varargin{vap_multirad}(1); else error(['CircularHough_Grd: ''multirad'' has to be ', ... 'within the range [0.1, 1]']); end end vap_fltr4accum = 4; % filter for smoothing the accumulation array if nargin > (1 + vap_fltr4accum), if isnumeric(varargin{vap_fltr4accum}) && ... ndims(varargin{vap_fltr4accum}) == 2 && ... all(size(varargin{vap_fltr4accum}) >= 3), fltr4accum = varargin{vap_fltr4accum}; else error(['CircularHough_Grd: ''fltr4accum'' has to be ', ... 'a 2-D matrix with a minimum size of 3-by-3']); end else % Default filter (5-by-5) fltr4accum = ones(5,5); fltr4accum(2:4,2:4) = 2; fltr4accum(3,3) = 6; end func_compu_cen = ( nargout > 1 ); func_compu_radii = ( nargout > 2 ); % Reserved parameters dbg_on = false; % debug information dbg_bfigno = 4; if nargout > 3, dbg_on = true; end %%%%%%%% Building accumulation array %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Convert the image to single if it is not of % class float (single or double) img_is_double = isa(img, 'double'); if ~(img_is_double || isa(img, 'single')), imgf = single(img); end % Compute the gradient and the magnitude of gradient if img_is_double, [grdx, grdy] = gradient(img); else [grdx, grdy] = gradient(imgf); end grdmag = sqrt(grdx.^2 + grdy.^2); % Get the linear indices, as well as the subscripts, of the pixels % whose gradient magnitudes are larger than the given threshold grdmasklin = find(grdmag > prm_grdthres); [grdmask_IdxI, grdmask_IdxJ] = ind2sub(size(grdmag), grdmasklin); %%%%%%%%%%%%%%%%%%is gradient_image SAME??????????%%%%%%%%%%%%%% [Height Width]=size(img); NEW_I=zeros(Height,Width); H=size(grdmask_IdxI); for i=1:H NEW_I(grdmask_IdxI(i,1),grdmask_IdxJ(i,1))=255; end imshow(NEW_I); % Compute the linear indices (as well as the subscripts) of % all the votings to the accumulation array. % The Matlab function 'accumarray' accepts only double variable, % so all indices are forced into double at this point. % A row in matrix 'lin2accum_aJ' contains the J indices (into the % accumulation array) of all the votings that are introduced by a % same pixel in the image. Similarly with matrix 'lin2accum_aI'. rr_4linaccum = double( prm_r_range ); linaccum_dr = [ (-rr_4linaccum(2) + 0.5) : -rr_4linaccum(1) , ... (rr_4linaccum(1) + 0.5) : rr_4linaccum(2) ]; %记录圆心的横坐标 lin2accum_aJ = floor( ... double(grdx(grdmasklin)./grdmag(grdmasklin)) * linaccum_dr + ... repmat( double(grdmask_IdxJ)+0.5 , [1,length(linaccum_dr)] ) ... ); %记录圆心的纵坐标 lin2accum_aI = floor( ... double(grdy(grdmasklin)./grdmag(grdmasklin)) * linaccum_dr + ... repmat( double(grdmask_IdxI)+0.5 , [1,length(linaccum_dr)] ) ... ); % Clip the votings that are out of the accumulation array %符合圆心条件的位置为1 mask_valid_aJaI = ... lin2accum_aJ > 0 & lin2accum_aJ < (size(grdmag,2) + 1) & ... lin2accum_aI > 0 & lin2accum_aI < (size(grdmag,1) + 1); mask_valid_aJaI_reverse = ~ mask_valid_aJaI; %记录圆心的实际横坐标 lin2accum_aJ = lin2accum_aJ .* mask_valid_aJaI + mask_valid_aJaI_reverse; %记录圆心的实际纵坐标 lin2accum_aI = lin2accum_aI .* mask_valid_aJaI + mask_valid_aJaI_reverse; clear mask_valid_aJaI_reverse; % Linear indices (of the votings) into the accumulation array lin2accum = sub2ind( size(grdmag), lin2accum_aI, lin2accum_aJ ); lin2accum_size = size( lin2accum ); lin2accum = reshape( lin2accum, [numel(lin2accum),1] ); %把圆心横纵坐标对应成原图的index % Weights of the votings, currently using the gradient maginitudes % but in fact any scheme can be used (application dependent) weight4accum = ... repmat( double(grdmag(grdmasklin)) , [lin2accum_size(2),1] ) .* ... mask_valid_aJaI(:); %weight4accum(:,1) clear mask_valid_aJaI; % Build the accumulation array using Matlab function 'accumarray' accum = accumarray( lin2accum , weight4accum ); accum = [ accum ; zeros( numel(grdmag) - numel(accum) , 1 ) ]; accum = reshape( accum, size(grdmag) ); %每个位置是圆心的投票 %%%%%%%% Locating local maxima in the accumulation array %%%%%%%%%%%% % Stop if no need to locate the center positions of circles if ~func_compu_cen, return; end clear lin2accum weight4accum; % Parameters to locate the local maxima in the accumulation array % -- Segmentation of 'accum' before locating LM prm_useaoi = true; prm_aoithres_s = 2; prm_aoiminsize = floor(min([ min(size(accum)) * 0.25, ... prm_r_range(2) * 1.5 ])); % -- Filter for searching for local maxima prm_fltrLM_s = 1.35; prm_fltrLM_r = ceil( prm_fltrLM_R * 0.6 ); prm_fltrLM_npix = max([ 6, ceil((prm_fltrLM_R/2)^1.8) ]); % -- Lower bound of the intensity of local maxima prm_LM_LoBndRa = 0.2; % minimum ratio of LM to the max of 'accum' % Smooth the accumulation array fltr4accum = fltr4accum / sum(fltr4accum(:)); accum = filter2( fltr4accum, accum ); %对圆心投票的结果进行了滤波 %imwrite(accum,'accum_TestImg_CHT_a3_filter.png','png'); % Select a number of Areas-Of-Interest from the accumulation array if prm_useaoi, % Threshold value for 'accum' prm_llm_thres1 = prm_grdthres * prm_aoithres_s; % Thresholding over the accumulation array accummask = ( accum > prm_llm_thres1 ); % Segmentation over the mask %imwrite(accummask,'accummask.png','png'); [accumlabel, accum_nRgn] = bwlabel( accummask, 8 ); %imwrite(accumlabel,'accumlabel.png','png'); % Select AOIs from segmented regions accumAOI = ones(0,4); for k = 1 : accum_nRgn, accumrgn_lin = find( accumlabel == k ); [accumrgn_IdxI, accumrgn_IdxJ] = ... ind2sub( size(accumlabel), accumrgn_lin ); rgn_top = min( accumrgn_IdxI ); rgn_bottom = max( accumrgn_IdxI ); rgn_left = min( accumrgn_IdxJ ); rgn_right = max( accumrgn_IdxJ ); % The AOIs selected must satisfy a minimum size if ( (rgn_right - rgn_left + 1) >= prm_aoiminsize && ... (rgn_bottom - rgn_top + 1) >= prm_aoiminsize ), accumAOI = [ accumAOI; ... rgn_top, rgn_bottom, rgn_left, rgn_right ]; end end else % Whole accumulation array as the one AOI accumAOI = [1, size(accum,1), 1, size(accum,2)]; end % Thresholding of 'accum' by a lower bound prm_LM_LoBnd = max(accum(:)) * prm_LM_LoBndRa; % Build the filter for searching for local maxima fltr4LM = zeros(2 * prm_fltrLM_R + 1); [mesh4fLM_x, mesh4fLM_y] = meshgrid(-prm_fltrLM_R : prm_fltrLM_R); mesh4fLM_r = sqrt( mesh4fLM_x.^2 + mesh4fLM_y.^2 ); fltr4LM_mask = ... ( mesh4fLM_r > prm_fltrLM_r & mesh4fLM_r <= prm_fltrLM_R ); fltr4LM = fltr4LM - ... fltr4LM_mask * (prm_fltrLM_s / sum(fltr4LM_mask(:))); if prm_fltrLM_R >= 4, fltr4LM_mask = ( mesh4fLM_r < (prm_fltrLM_r - 1) ); else fltr4LM_mask = ( mesh4fLM_r < prm_fltrLM_r ); end fltr4LM = fltr4LM + fltr4LM_mask / sum(fltr4LM_mask(:)); % **** Debug code (begin) if dbg_on, dbg_LMmask = zeros(size(accum)); end % **** Debug code (end) % For each of the AOIs selected, locate the local maxima circen = zeros(0,2); for k = 1 : size(accumAOI, 1), aoi = accumAOI(k,:); % just for referencing convenience % Thresholding of 'accum' by a lower bound accumaoi_LBMask = ... ( accum(aoi(1):aoi(2), aoi(3):aoi(4)) > prm_LM_LoBnd ); % Apply the local maxima filter candLM = conv2( accum(aoi(1):aoi(2), aoi(3):aoi(4)) , ... fltr4LM , 'same' ); candLM_mask = ( candLM > 0 ); % Clear the margins of 'candLM_mask' candLM_mask([1:prm_fltrLM_R, (end-prm_fltrLM_R+1):end], :) = 0; candLM_mask(:, [1:prm_fltrLM_R, (end-prm_fltrLM_R+1):end]) = 0; % **** Debug code (begin) if dbg_on, dbg_LMmask(aoi(1):aoi(2), aoi(3):aoi(4)) = ... dbg_LMmask(aoi(1):aoi(2), aoi(3):aoi(4)) + ... accumaoi_LBMask + 2 * candLM_mask; end % **** Debug code (end) % Group the local maxima candidates by adjacency, compute the % centroid position for each group and take that as the center % of one circle detected [candLM_label, candLM_nRgn] = bwlabel( candLM_mask, 8 ); for ilabel = 1 : candLM_nRgn, % Indices (to current AOI) of the pixels in the group candgrp_masklin = find( candLM_label == ilabel ); [candgrp_IdxI, candgrp_IdxJ] = ... ind2sub( size(candLM_label) , candgrp_masklin ); % Indices (to 'accum') of the pixels in the4 group candgrp_IdxI = candgrp_IdxI + ( aoi(1) - 1 ); candgrp_IdxJ = candgrp_IdxJ + ( aoi(3) - 1 ); candgrp_idx2acm = ... sub2ind( size(accum) , candgrp_IdxI , candgrp_IdxJ ); % Minimum number of qulified pixels in the group if sum(accumaoi_LBMask(candgrp_masklin)) < prm_fltrLM_npix, continue; end % Compute the centroid position candgrp_acmsum = sum( accum(candgrp_idx2acm) ); cc_x = sum( candgrp_IdxJ .* accum(candgrp_idx2acm) ) / ... candgrp_acmsum; cc_y = sum( candgrp_IdxI .* accum(candgrp_idx2acm) ) / ... candgrp_acmsum; circen = [circen; cc_x, cc_y]; end end % **** Debug code (begin) if dbg_on, figure(dbg_bfigno); imagesc(dbg_LMmask); axis image; title('Generated map of local maxima'); if size(accumAOI, 1) == 1, figure(dbg_bfigno+1); surf(candLM, 'EdgeColor', 'none'); axis ij; title('Accumulation array after local maximum filtering'); end end % **** Debug code (end) %%%%%%%% Estimation of the Radii of Circles %%%%%%%%%%%% % Stop if no need to estimate the radii of circles if ~func_compu_radii, varargout{1} = circen; return; end % Parameters for the estimation of the radii of circles fltr4SgnCv = [2 1 1]; fltr4SgnCv = fltr4SgnCv / sum(fltr4SgnCv); % Find circle's radius using its signature curve cirrad = zeros( size(circen,1), 1 ); for k = 1 : size(circen,1), % Neighborhood region of the circle for building the sgn. curve circen_round = round( circen(k,:) ); SCvR_I0 = circen_round(2) - prm_r_range(2) - 1; if SCvR_I0 < 1, SCvR_I0 = 1; end SCvR_I1 = circen_round(2) + prm_r_range(2) + 1; if SCvR_I1 > size(grdx,1), SCvR_I1 = size(grdx,1); end SCvR_J0 = circen_round(1) - prm_r_range(2) - 1; if SCvR_J0 < 1, SCvR_J0 = 1; end SCvR_J1 = circen_round(1) + prm_r_range(2) + 1; if SCvR_J1 > size(grdx,2), SCvR_J1 = size(grdx,2); end % Build the sgn. curve SgnCvMat_dx = repmat( (SCvR_J0:SCvR_J1) - circen(k,1) , ... [SCvR_I1 - SCvR_I0 + 1 , 1] ); SgnCvMat_dy = repmat( (SCvR_I0:SCvR_I1)' - circen(k,2) , ... [1 , SCvR_J1 - SCvR_J0 + 1] ); SgnCvMat_r = sqrt( SgnCvMat_dx .^2 + SgnCvMat_dy .^2 ); SgnCvMat_rp1 = round(SgnCvMat_r) + 1; f4SgnCv = abs( ... double(grdx(SCvR_I0:SCvR_I1, SCvR_J0:SCvR_J1)) .* SgnCvMat_dx + ... double(grdy(SCvR_I0:SCvR_I1, SCvR_J0:SCvR_J1)) .* SgnCvMat_dy ... ) ./ SgnCvMat_r; SgnCv = accumarray( SgnCvMat_rp1(:) , f4SgnCv(:) ); SgnCv_Cnt = accumarray( SgnCvMat_rp1(:) , ones(numel(f4SgnCv),1) ); SgnCv_Cnt = SgnCv_Cnt + (SgnCv_Cnt == 0); SgnCv = SgnCv ./ SgnCv_Cnt; % Suppress the undesired entries in the sgn. curve % -- Radii that correspond to short arcs SgnCv = SgnCv .* ( SgnCv_Cnt >= (pi/4 * [0:(numel(SgnCv_Cnt)-1)]') ); % -- Radii that are out of the given range SgnCv( 1 : (round(prm_r_range(1))+1) ) = 0; SgnCv( (round(prm_r_range(2))+1) : end ) = 0; % Get rid of the zero radius entry in the array SgnCv = SgnCv(2:end); % Smooth the sgn. curve SgnCv = filtfilt( fltr4SgnCv , [1] , SgnCv ); % Get the maximum value in the sgn. curve SgnCv_max = max(SgnCv); if SgnCv_max <= 0, cirrad(k) = 0; continue; end % Find the local maxima in sgn. curve by 1st order derivatives % -- Mark the ascending edges in the sgn. curve as 1s and % -- descending edges as 0s SgnCv_AscEdg = ( SgnCv(2:end) - SgnCv(1:(end-1)) ) > 0; % -- Mark the transition (ascending to descending) regions SgnCv_LMmask = [ 0; 0; SgnCv_AscEdg(1:(end-2)) ] & (~SgnCv_AscEdg); SgnCv_LMmask = SgnCv_LMmask & [ SgnCv_LMmask(2:end) ; 0 ]; % Incorporate the minimum value requirement SgnCv_LMmask = SgnCv_LMmask & ... ( SgnCv(1:(end-1)) >= (prm_multirad * SgnCv_max) ); % Get the positions of the peaks SgnCv_LMPos = sort( find(SgnCv_LMmask) ); % Save the detected radii if isempty(SgnCv_LMPos), cirrad(k) = 0; else cirrad(k) = SgnCv_LMPos(end); for i_radii = (length(SgnCv_LMPos) - 1) : -1 : 1, circen = [ circen; circen(k,:) ]; cirrad = [ cirrad; SgnCv_LMPos(i_radii) ]; end end end % Output varargout{1} = circen; varargout{2} = cirrad; if nargout > 3, varargout{3} = dbg_LMmask; end