/matlab/src/CircularHough_Grd.m
MATLAB | 651 lines | 280 code | 65 blank | 306 comment | 35 complexity | e0400ede69347d234274267e0bb26cde MD5 | raw file
- %> @file CircularHough_Grd.m
- %> @brief HOUGH??????
- %=============================================
- %> ????,??????????????
- %>
- 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 );
-
- parfor 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
-