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);
% 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
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] );
clear lin2accum_aI lin2accum_aJ;
% 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(:);
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 );
% 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
[accumlabel, accum_nRgn] = bwlabel( accummask, 8 );
% 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 the 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
