1## Copyright (C) 2017 Hartmut Gimpel <hg_code@gmx.de> 2## Copyright (C) 2018 Carnë Draug <carandraug@octave.org> 3## 4## This program is free software; you can redistribute it and/or modify 5## it under the terms of the GNU General Public License as published by 6## the Free Software Foundation; either version 3 of the License, or 7## (at your option) any later version. 8## 9## This program is distributed in the hope that it will be useful, 10## but WITHOUT ANY WARRANTY; without even the implied warranty of 11## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12## GNU General Public License for more details. 13## 14## You should have received a copy of the GNU General Public License 15## along with this program. If not, see <http://www.gnu.org/licenses/>. 16 17## -*- texinfo -*- 18## @deftypefn {Function File} {@var{centers} =} imfindcircles (@var{im}, @var{radius}) 19## @deftypefnx {Function File} {@var{centers} =} imfindcircles (@var{im}, @var{RadiusRange}) 20## @deftypefnx {Function File} {@var{centers} =} imfindcircles (@dots{}, @var{property}, @var{value}) 21## @deftypefnx {Function File} {[@var{centers}, @var{radii}, @var{strengths}] =} imfindcircles (@dots{}) 22## Find circles in image using the circular Hough transform. 23## 24## This function finds circles in a given 2d image @var{im}. 25## The image data can be grayscale, rgb color, or binary. 26## The search is done only for circular shapes of approximatly 27## the given @var{radius} (single value) or @var{RadiusRange} 28## (a two element vector [r_min, r_max]). 29## 30## The output @var{centers} is a two column matrix, 31## where each row contains the (fractional) x and y coordinates of 32## the center of one found circle. The output rows are ordered 33## according the @var{strengths} of the circles, the strongest 34## circle first. 35## 36## The output @var{radii} is a column vector of the (fractional) 37## radii of the found circles, in the same order as the circle centers. 38## If @var{radius} is specified, instead of @var{RadiusRange}, then 39## all values in @var{radii} are the same as @var{radius}. 40## 41## The output @var{strengths} is a measure of how "strong" 42## the corresponding circle is. (A sharp, smooth and full circle gives 43## roughly a value of one. The minimum value is zero.) 44## 45## Additionally the following optional property-value-pairs can be used: 46## @table @var 47## @item ObjectPolarity 48## Either "bright" circles on darker background can be searched (default), 49## or "dark" circles on brighter background. 50## 51## @item Method 52## The default method "PhaseCode" is faster as the 53## (currently unimplemented) method "TwoStage". 54## For literature and details on the algorithm, see comments 55## in the code. 56## 57## @item Sensitivity 58## A value between 0 and 1 to say how strong 59## a circle must be in order to be found. With a value 60## of 1 no circles are discarded, default value is 0.85. 61## (Use bigger values of Sensitivity if your circles are 62## not properly found, until unwanted "ghost" circles 63## start to appear.) 64## 65## @item EdgeThreshold 66## A value between 0 and 1 to say how strong 67## an edge point of the image must be, to be 68## considered as a possible point on a circle circumference. 69## With a value of 0 all possible 70## edge points are considered, with a value of 1 only the 71## edge point with the stronges gradient will be considered. 72## As default value the output of @code{graythresh} (Otsu's threshold) 73## of the gradient image is taken. 74## @end table 75## 76## Notes: 77## @itemize @bullet 78## @item 79## Only center points inside the image region can be found. 80## @item 81## For concentric cirles the output is unpredictable. 82## @item 83## For optimal speed, keep the radius range and the image size 84## as small as possible. 85## @item 86## For radius values above 100 pixels the sensitivity and accuracy 87## of this algorithm starts to decrease. 88## @item 89## For big radius ranges the sensitivity decreases. Try to keep 90## r_max < 3 * r_min . 91## @item 92## RGB color images will automatically be converted to grayscale 93## using the @code{rgb2gray} function. 94## @item 95## Binary images will automatically be converted to grayscale 96## and be slightly smoothed before processing. 97## @end itemize 98## 99## Compatibility note: The @var{centers} and @var{radii} outputs 100## have good compatibility to the Matlab outputs. The @var{strengths} 101## outputs differ a bit, but have mostly the same ordering. 102## Currently only the (default) Matlab algorithm "PhaseCode" is implemented 103## in Octave. 104## 105## @seealso{hough, hough_circle} 106## @end deftypefn 107 108## Algorithm: 109## The following papers and books were used for the 110## implementation of this function: 111## 112## [1] (for general information on the circle Hough transform 113## and the basic algorithmic approach) 114## E. R. Davies: "Computer & Machine Vision" 115## Academic Press (2012), 4th edition 116## chapter 12 "Circle and Ellipse Detection" 117## [2] (for the 'phase code' algorithm to search for a radius range 118## with just a 2-dimensional accumulator array) 119## T. J. Artherton and D. J. Kerbyson: 120## "Size invariant circle detection" 121## Image and Vision Computing 17 (1999), p. 795-803. 122## [3] (for 'state of the art' peak finding in the circular Hough 123## transform accumulator array) 124## C. Zhang, F. Huber, M. Knop, F. A. Hamprecht: 125## "Yest cell detection and segmentation in bright field microscopy" 126## IEEE Int. Symp. Biomed. Imag. (ISBI), April 2014, p. 1267-1270. 127## 128## A more detailed description of the individual steps of the 129## algorithm is given in the following code itself. 130 131function [centers, radii, strengths] = imfindcircles (im, radiusrange, varargin) 132 133 if (nargin < 2 || nargin > 10 || mod (numel (varargin), 2) != 0) 134 print_usage (); 135 endif 136 137 p = inputParser (); 138 p.FunctionName = "imfindcircles"; 139 p.addParamValue ("ObjectPolarity", "bright", 140 @(x) any (strcmpi (x, {"bright", "dark"}))); 141 p.addParamValue ("Method", "PhaseCode", 142 @(x) any (strcmpi (x, {"PhaseCode", "TwoStage"}))); 143 p.addParamValue ("Sensitivity", 0.85, 144 @(x) isnumeric (x) && isreal (x) && isscalar (x)); 145 146 ## The default value of EdgeThreshold must be computed later. 147 p.addParamValue ("EdgeThreshold", [], 148 @(x) isnumeric (x) && isreal (x) && isscalar (x)); 149 p.parse (varargin{:}); 150 151 dark_circles = strcmpi (p.Results.ObjectPolarity, "dark"); 152 sensitivity = p.Results.Sensitivity; 153 edge_thresh = p.Results.EdgeThreshold; 154 155 if (! isimage (im) || ! any (ndims (im) == [2, 3])) 156 error ("imfindcircles: IM must be a logical or numeric 2d or 3d array"); 157 elseif (ndims (im) == 3 && size (im, 3) != 3) 158 error ("imfindcircles: the 3d image IM must be a RGB image"); 159 endif 160 161 if (all (numel (radiusrange) != [1 2])) 162 error ("imfindcircles: RADIUS or RADIUSRANGE must be a vector of length 1 or 2"); 163 elseif (! isnumeric (radiusrange) || any (radiusrange <= 0)) 164 error ("imfindcircles: RADIUS or RADIUSRANGE must be positive") 165 elseif (numel (radiusrange) == 2 && radiusrange(1) > radiusrange(2)) 166 error ("imfindcircles: RADIUSRANGE(1) must be smaller than RADIUSRANGE(2)") 167 endif 168 169 if (strcmpi (p.Results.Method, "TwoStage")) 170 error ("imfindcirles: the 'TwoStage' method is not yet implemented. Try 'PhaseCode' instead."); 171 endif 172 173 if (sensitivity < 0 || sensitivity > 1) 174 error ("imfindcircles: 'Sensitivity' must be between 0 and 1"); 175 endif 176 177 if (! isempty (edge_thresh) && (edge_thresh < 0 || edge_thresh > 1)) 178 error ("imfindcircles: 'EdgeThreshold' must be between 0 and 1"); 179 endif 180 181 H = cht_accumulator (im, radiusrange, edge_thresh, dark_circles); 182 [centers, centers_ind, strengths] = cht_centers (H, sensitivity); 183 184 ## Output values must be sorted by strength. 185 [strengths, idx_sorted] = sort (strengths, "descend"); 186 centers = centers(idx_sorted,:); 187 188 if (isargout (2)) 189 if (isscalar (radiusrange)) 190 radii = repmat (radiusrange, numel (strengths), 1); 191 else 192 radii = cht_radii (H(centers_ind), radiusrange); 193 ## Output values must be sorted by strength. 194 radii = radii(idx_sorted); 195 endif 196 endif 197endfunction 198 199 200## Convert image from RGB or logical types into a floating point 2D 201## image that can be used for CHT. 202function im = cht_image_preparation (im) 203 if (isinteger (im)) 204 im = im2double (im); 205 endif 206 207 if (size (im, 3) == 3) 208 im = rgb2gray (im); 209 endif 210 211 ## The Matlab help page mentions "preprocessing" of 212 ## binary input images to "improve the accuracy". 213 ## Unsure what they do here. 214 ## We'll just do a little (carefully rotation symmetrical) 215 ## smoothing as grayscale images, because we need 216 ## the accurate gradient directions later on. (This wouldn't 217 ## work properly with the binary images itself.) 218 if (islogical (im)) 219 im = im2single (im); 220 gauss_filter = fspecial ("gaussian", [9 9], 1.5); 221 im = imfilter (im, gauss_filter, "replicate"); 222 endif 223endfunction 224 225 226## calculate the intensity gradients values at edge pixels 227## using the Sobel operator (see ref. [1], chapter 12.2) 228## (we'll always call the first coordinate "x" here) 229function [G0, GxE, GyE, E_ind] = edge_intensity_gradients (im, edge_thresh) 230 sobel_x = fspecial ("sobel"); 231 sobel_y = sobel_x'; 232 Gx = imfilter (im, sobel_x, "replicate"); 233 Gy = imfilter (im, sobel_y, "replicate"); 234 G = sqrt (Gx.^ 2 + Gy.^ 2); # sqrt is expensive, but yields higher angle precision 235 236 ## get rid of empty gradient images, they would cause trouble later on: 237 Gmax = max (G(:)); 238 if (Gmax == 0) 239 E_ind = [](:); 240 else 241 ## Find the edge pixels by thresholding the gradient image. If 242 ## not set, estimate threshold with Otsu's algorithm. 243 if (isempty (edge_thresh)) 244 edge_thresh = graythresh (G ./ Gmax); 245 endif 246 E_ind = find (G > (edge_thresh * Gmax)); # edge pixels indices 247 endif 248 G0 = G(E_ind); 249 GxE = Gx(E_ind); 250 GyE = Gy(E_ind); 251endfunction 252 253 254function [xs, ys] = candidate_centers (im, R, edge_thresh) 255 [G0, GxE, GyE, E_ind] = edge_intensity_gradients (im, edge_thresh); 256 [Ex, Ey] = ind2sub (size (im), E_ind); 257 258 ## generate stroke pixels (xs, ys): 259 ## This does not involve sin and cos calculations, 260 ## which helps for speed (see ref. [1].) 261 ## (automatic broadcasting here, 262 ## R is a row vector, the G's are column vectors.) 263 xs = Ex - R .* (GxE ./ G0); # eq. 12.3 (ref. [1]) 264 ys = Ey - R .* (GyE ./ G0); # eq. 12.4 (ref. [1]) 265 ## round to integer pixel positions: 266 xs = round (xs); 267 ys = round (ys); 268endfunction 269 270 271## visit all edge pixels and generate 272## a "spoke" in the circular Hough transform: 273## A spoke is a line of pixels which are a distance r1 to r2 274## away from the edge point, this distance is taken in 275## the edge direction (gradient) of this very edge point. 276## The resulting points on the spokes are candidates for circle 277## center points. And each gradient pixel can vote for 278## one center pixel candidate (for each radius value). 279## All spoke points are summed up in the accumulator array H. 280## (see equations 12.3 and 12.4 for xs and ys in ref. [1]) 281function H = cht_accumulator (im, radiusrange, edge_thresh, dark_circles) 282 im = cht_image_preparation (im); 283 284 if (isscalar (radiusrange)) 285 R = radiusrange; 286 else 287 ## Range of circle radii. Step size of 0.5 covers all integer 288 ## circle diameters. 289 R = radiusrange(1):0.5:radiusrange(2); 290 endif 291 292 ## code the circle polarity in the radius sign: 293 ## (different direction of "spoke", see ref. [1], chapter 12.2) 294 if (dark_circles) 295 [xs, ys] = candidate_centers (im, -R, edge_thresh); 296 else 297 [xs, ys] = candidate_centers (im, R, edge_thresh); 298 endif 299 ns = rows (xs); 300 301 ## limit pixel position and code values to the image size: 302 xs = xs(:); 303 ys = ys(:); 304 idx_outside = (xs < 1 | xs > rows (im) | ys < 1 | ys > columns (im)); 305 xs(idx_outside) = []; 306 ys(idx_outside) = []; 307 308 if (isscalar (radiusrange)) 309 codes = 1 ./ R; # phase coding is not useful for a single radius value 310 else 311 ## pre-calculate the phases and (complex) code values for all radii: 312 ## We use the "log phase coding", which is best, according to 313 ## ref. [2]. 314 logR = log (R); 315 phase = 2 .* pi .* ((logR - logR(1)) 316 ./ (logR(end) - logR(1))); # eq. 8 (ref. [2]) 317 code = exp (i .* phase); # eq. 5 (ref. [2]) 318 code ./= R; # eq. 11 (ref. [2]) , radius normalization 319 ## -> a full and perfect circle will give a value of 2pi in H 320 321 codes = repmat (code, ns, 1)(:); 322 codes(idx_outside) = []; 323 endif 324 325 ## add the code value of all those spoke pixels 326 ## to the corresponding pixel position in H: 327 H = accumarray ([xs, ys], codes, size (im)); 328endfunction 329 330 331## Find local (regional) maximas in H, because they are the circle 332## center pixels. This part of the algorithm is taken from ref. [3], 333## section 2.2 "Detecting cell center candidates". 334function [centers, centers_ind, strengths] = cht_centers (H, sensitivity) 335 ## Take the absolute value of the H accumulator array which may be 336 ## complex if we had a radius range and phase coded annulus. 337 H = abs (H); 338 339 ## do some smoothing before searching for peak positions 340 ## use a (rotation symmetric) gaussian filter for this (see ref. [3]) 341 gauss_filter = fspecial ("gaussian", [9 9], 2); # a bigger filter might join nearby centers 342 H = imfilter (H, gauss_filter, 0); # introduces small shifts of centers near the edge, 343 # but padding "replicate" would be worse 344 345 ## define threshold to suppress smaller peaks in H: 346 peak_thresh = 1 - sensitivity; 347 peak_thresh *= 2*pi; # because a full circle can add up to 2pi in H 348 349 ## use the h-maxima transform to suppress maximas with a too low height: 350 Hbig = imhmax (H, peak_thresh); 351 ## find the regional maxima of those smoothed big peaks in H: 352 Hmax_BW = imregionalmax (Hbig ./ max (Hbig(:))); 353 354 ## Use those maxima as well as their surrounding peak 355 ## to calculate a weighted average for the peak position: 356 ## (This gives the circle center positions with the necessary 357 ## subpixel accuracy.) 358 props = regionprops (Hmax_BW, Hbig, "WeightedCentroid"); 359 centers = cell2mat ({props.WeightedCentroid}'); 360 361 ## for the "strength" return value, take the (smoothed) H: 362 ## (Matlab seems to do this a bit different.) 363 centers_ind = sub2ind (size (H), round (centers(:,2)), round (centers(:,1))); 364 strengths = H(centers_ind) ./ (2*pi); # normalize a "full" circle to 1 365endfunction 366 367 368## calculate the circle radius from the complex phase in H: 369## (i.e. undo the phase coding) 370function radii = cht_radii (H, radiusrange) 371 r1 = radiusrange(1); 372 r2 = radiusrange(2); 373 c_phase = arg (H); 374 c_phase(c_phase < 0) += (2*pi); 375 radii = exp ((c_phase ./ (2*pi) .* (log (r2) - log (r1))) + log (r1)); 376endfunction 377 378 379%!shared im0, rgb0, im1 380%! im0 = [0 0 0 0 0; 381%! 0 1 2 1 0; 382%! 0 2 5 2 0; 383%! 0 1 2 1 0; 384%! 0 0 0 0 0]; 385%! rgb0 = cat (3, im0, 3.*im0, 2.*im0); 386%! im1 = zeros (20); 387%! im1(2:6, 5:9) = 1; 388%! im1(13:19, 13:19) = 1; 389 390%!function image = circlesimage (numx, numy, centersx, centersy, rs, values) 391%! ## create an image with circles of given parameters 392%! num = length (centersx); 393%! image = zeros (numy, numx); 394%! [indy, indx] = meshgrid (1:numx, 1:numy); 395%! for n = 1:num 396%! centerx = centersx(n); 397%! centery = centersy(n); 398%! r = rs(n); 399%! value = values(n); 400%! dist_squared = (indx - centerx).^ 2 + (indy - centery).^ 2; 401%! image(dist_squared <= (r-0.5)^2) = value; 402%! endfor 403%!endfunction 404 405 406## test input syntax: 407%!error imfindcircles () 408%!error imfindcircles (im0) 409%!error imfindcircles (im0, [1 2 3]) 410%!error imfindcircles (im0, -3) 411%!error imfindcircles (im0, 4+2*i) 412%!error imfindcircles (ones (5,5,4), 2) 413%!error imfindcircles (ones (5,5,5,5), 2) 414%!error imfindcircles (im0, [2 1]) 415%!error imfindcircles (im0, 2, "rubbish") 416%!error imfindcircles (im0, 2, "more", "rubbish") 417%!error imfindcircles (im0, 2, "ObjectPolarity", "rubbish") 418%!error imfindcircles (im0, 2, "ObjectPolarity", 5) 419%!error imfindcircles (im0, 2, "ObjectPolarity") 420%!error imfindcircles (im0, 2, "Method", "rubbish") 421%!error imfindcircles (im0, 2, "Method", 5) 422%!error imfindcircles (im0, 2, "Method") 423%!error imfindcircles (im0, 2, "Sensitivity", "rubbish") 424%!error imfindcircles (im0, 2, "Sensitivity") 425%!error imfindcircles (im0, 2, "Sensitivity", -0.1) 426%!error imfindcircles (im0, 2, "Sensitivity", 1.1) 427%!error imfindcircles (im0, 2, "Sensitivity", [0.1 0.2]) 428%!error imfindcircles (im0, 2, "EdgeThreshold", "rubbish") 429%!error imfindcircles (im0, 2, "EdgeThreshold") 430%!error imfindcircles (im0, 2, "EdgeThreshold", -0.1) 431%!error imfindcircles (im0, 2, "EdgeThreshold", 1.1) 432%!error imfindcircles (im0, 2, "EdgeThreshold", [0.1 0.2]) 433%!error imfindcircles (im0, 2, "EdgeThreshold", 0.1, "ObjectPolarity", "bright", 434%! "Sensitivity", 0.3, "Method", "PhaseCode", "more", 1) 435 436%!test # none of this should fail 437%! imfindcircles (im0, 2); 438%! imfindcircles (im0, [1 2]); 439%! imfindcircles (logical (im0), 2); 440%! imfindcircles (logical (im0), [1 2]); 441%! imfindcircles (rgb0, 2); 442%! imfindcircles (rgb0, [1 2]); 443%! imfindcircles (uint8 (im0), 2); 444%! imfindcircles (uint8 (im0), [1 2]); 445%! imfindcircles (im0, 2, "ObjectPolarity", "bright"); 446%! imfindcircles (im0, 2, "ObjectPolarity", "dark"); 447%! imfindcircles (im0, 2, "Method", "PhaseCode"); 448%! imfindcircles (im0, 2, "Sensitivity", 0.5); 449%! imfindcircles (im0, 2, "EdgeThreshold", 0.5); 450%! imfindcircles (im0, 2, "ObjectPolarity", "bright", "Method", "PhaseCode"); 451%! imfindcircles (im0, 2, "ObjectPolarity", "bright", "Sensitivity", 0.3, 452%! "Method", "PhaseCode"); 453%! imfindcircles (im0, 2, "EdgeThreshold", 0.1, "ObjectPolarity", "bright", 454%! "Sensitivity", 0.3, "Method", "PhaseCode"); 455 456 457## output class, number and shape: 458%!test 459%! centers = imfindcircles (im1, 2); 460%! assert (size (centers, 2), 2) 461%! assert (class (centers), "double") 462 463%!test 464%! [centers, radii] = imfindcircles (im1, [1 5]); 465%! assert (size (centers, 2), 2) 466%! assert (size (radii, 2), 1) 467%! assert (class (radii), "double") 468 469%!test 470%! [centers, radii, strengths] = imfindcircles (im1, [1 5]); 471%! assert (size (strengths, 2), 1) 472%! assert (class (strengths), "double") 473 474%!error [a b c d] = imfindcircles (im0, 2); 475 476 477## test calculation results: 478%!test ## sub-pixel accuracy of circle center 479%! xs = [95.7]; 480%! ys = [101.1]; 481%! rs = [50]; 482%! vals = [0.5]; 483%! im = circlesimage (200, 200, xs, ys, rs, vals); 484%! filt = ones (3) ./ 9; 485%! im = imfilter (im, filt); 486%! [centers, radii] = imfindcircles (im, [40 60]); 487%! assert (centers, [101.1, 95.7], 0.1); 488%! assert (radii, 50, 1); 489 490%!test 491%! ## specificity to circular shapes and strengths output value 492%! xs = [100 202]; 493%! ys = [101, 203]; 494%! rs = [40, 41]; 495%! vals = [0.8, 0.9]; 496%! im = circlesimage (300, 300, xs, ys, rs, vals); 497%! filt = ones (3) ./ 9; 498%! im = imfilter (im, filt); 499%! im(30:170, 50:100) = 0; 500%! im(20:120, 180:280) = 1; 501%! [centers, radii, strengths] = imfindcircles (im, [30 50], "Sensitivity", 0.9); 502%! assert (size (centers), [2 2]); 503%! assert (centers, [203, 202; 101, 100], 1.5); 504%! assert (radii, [40; 41], 2.5); 505%! assert (strengths(1) / strengths(2) > 1.8, true); 506 507%!test # radius range parameter & dark circles 508%! xs = [50, 420, 180]; 509%! ys = [80, 100, 200]; 510%! rs = [35, 30, 40]; 511%! vals = [0.7, 0.8, 0.9]; 512%! im = circlesimage (300, 500, xs, ys, rs, vals); 513%! filt = ones (3) ./ 9; 514%! im = imfilter (im, filt); 515%! [centers1, radii1] = imfindcircles (im, [28 36]); 516%! [centers2, radii2] = imfindcircles (im, [28 42]); 517%! assert (size (centers1), [2 2]); 518%! assert (centers1, [100 420; 80 50], 0.2); 519%! assert (radii1, [30; 35], 2); 520%! assert (size (centers2), [3 2]); 521%! im_dark = 1-im; 522%! [centers_dark, radii_dark, strengths_dark] = imfindcircles (im_dark, [25 42], "ObjectPolarity", "dark"); 523%! assert (sortrows (centers_dark), [80 50; 100 420; 200 180], 0.2); 524%! assert (sortrows (radii_dark), [30; 35; 40], 1); 525 526%!test # ability to find circles with big radius 527%! xs = [111, 555, 341]; 528%! ys = [222, 401, 161]; 529%! rs = [45, 50, 150]; 530%! vals = [0.6, 0.8, 0.7]; 531%! im = circlesimage (400, 701, xs, ys, rs, vals); 532%! [centers, radii] = imfindcircles (im, [140 160], "Sensitivity", 0.98); 533%! assert (centers, [161, 341], 0.2); 534%! assert (radii, 150, 1); 535 536%!test # overlapping circles 537%! xs = [105, 155]; 538%! ys = [202, 221]; 539%! rs = [45, 50]; 540%! vals = [0.5, 0.8]; 541%! im = circlesimage(385, 422, xs, ys, rs, vals); 542%! filt = ones (3) ./ 9; 543%! im = imfilter (im, filt); 544%! [centers, radii] = imfindcircles (im, [30 80]); 545%! assert (centers, [221, 155; 202, 105], 0.5); 546%! assert (radii, [50; 45], 1); 547 548%!test # overlapping circles, only 10 pixels apart 549%! xs = [155, 155]; 550%! ys = [175, 157]; 551%! rs = [50, 50]; 552%! vals = [0.7, 0.8]; 553%! im = circlesimage (300, 300, xs, ys, rs, vals); 554%! filt = ones (3) ./ 9; 555%! im = imfilter (im, filt); 556%! [centers, radii] = imfindcircles (im, [30 80], "Sensitivity", 0.95); 557%! assert (centers, [157, 155; 175, 155], 1); 558%! assert (radii, [50; 50], 1); 559 560%!test # edge threshold parameter 561%! xs = [100 202]; 562%! ys = [101, 203]; 563%! rs = [40, 41]; 564%! vals = [0.1, 0.9]; 565%! im = circlesimage (300, 300, xs, ys, rs, vals); 566%! filt = ones (3) ./ 9; 567%! im= imfilter (im, filt); 568%! [centers_auto, radii_auto] = imfindcircles (im, [30 50]); 569%! [centers_0, radii_0] = imfindcircles (im, [30 50], "EdgeThreshold", 0); 570%! [centers_05, radii_05] = imfindcircles (im, [30 50], "EdgeThreshold", 0.5); 571%! assert (centers_auto, [203, 202], 0.2); 572%! assert (radii_auto, 41, 1); 573%! assert (centers_0, [101, 100; 203, 202], 0.2); 574%! assert (radii_0, [40; 41], 1); 575%! assert (centers_05, [203, 202], 0.2); 576%! assert (radii_05, 41, 1); 577 578%!demo 579%! ## First generate an input image: 580%! model = [ 1.0 0.2 0.2 0.2 0.5 0 581%! 1.0 0.3 0.3 -0.1 -0.2 0 582%! -0.5 0.7 0.7 -0.5 0.5 0]; 583%! im = phantom (model); 584%! im(170:230,170:230) = 1; 585%! im = imfilter (im, fspecial ("average", 3)); 586%! im = imnoise (im, "salt & pepper"); 587%! imshow (im); 588%! 589%! ## Find and show circles with radius between 20 and 50: 590%! [centers, radii] = imfindcircles (im, [20 50]); 591%! viscircles (centers, radii) 592%! title ("found circles in red") 593