octave_packages/image-1.0.15/hough_circle.m

   1 ## Copyright (C) 2008 Soren Hauberg
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 3 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; if not, write to the Free Software
  15 ## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  16
  17 ## -*- texinfo -*-
  18 ## @deftypefn {Function File} @var{accum}= hough_circle (@var{bw}, @var{r})
  19 ## Perform the Hough transform for circles with radius @var{r} on the
  20 ## black-and-white image @var{bw}.
  21 ##
  22 ## As an example, the following shows how to compute the Hough transform for circles
  23 ## with radius 3 or 7 in the image @var{im}
  24 ## @example
  25 ## bw = edge(im);
  26 ## accum = hough_circle(bw, [3, 7]);
  27 ## @end example
  28 ## If @var{im} is an NxM image @var{accum} will be an NxMx2 array, where
  29 ## @var{accum}(:,:,1) will contain the Hough transform for circles with
  30 ## radius 3, and @var{accum}(:,:,2) for radius 7. To find good circles you
  31 ## now need to find local maximas in @var{accum}, which can be a hard problem.
  32 ## If you find a local maxima in @var{accum}(row, col, 1) it means that a
  33 ## good circle exists with center (row,col) and radius 3.
  34 ##
  35 ## @seealso{houghtf}
  36 ## @end deftypefn
  37
  38 function accum = hough_circle(bw, r)
  39   ## Check input arguments
  40   if (nargin != 2)
  41     error("hough_circle: wrong number of input arguments");
  42   endif
  43
  44   if (!ismatrix(bw) || ndims(bw) != 2)
  45     error("hough_circle: first arguments must be a 2-dimensional matrix");
  46   endif
  47
  48   if (!isvector(r) || !isreal(r) || any(r<0))
  49     error("hough_circle: radius arguments must be a positive vector or scalar");
  50   endif
  51
  52   ## Create the accumulator array.
  53   accum = zeros(size(bw,1), size(bw,2), length(r));
  54
  55   ## Find the pixels we need to look at
  56   [R, C] = find(bw);
  57
  58   ## Iterate over different radius
  59   for j = 1:length(r)
  60     rad = r(j);
  61
  62     ## Compute a filter containing the circle we're looking for.
  63     circ = circle(rad);
  64
  65     ## Iterate over all interesting image points
  66     for i =1:length(R)
  67       row = R(i);
  68       col = C(i);
  69
  70       ## Compute indices for the accumulator array
  71       a_rows = max(row-rad,1) : min(row+rad, size(accum,1));
  72       a_cols = max(col-rad,1) : min(col+rad, size(accum,2));
  73
  74       ## Compute indices for the circle array (the filter)
  75       c_rows = max(rad-row+2,1) : min(rad-row+1+size(accum,1), size(circ,1));
  76       c_cols = max(rad-col+2,1) : min(rad-col+1+size(accum,2), size(circ,2));
  77
  78       ## Update the accumulator array
  79       accum( a_rows, a_cols, j ) += circ ( c_rows, c_cols );
  80     endfor
  81   endfor
  82 endfunction
  83
  84 ## Small auxilary function that creates an (2r+1)x(2r+1) image containing
  85 ## a circle with radius r and center (r+1, r+1).
  86 function circ = circle(r)
  87   circ = zeros(round(2*r+1));
  88   col = 1:size(circ,2);
  89   for row=1:size(circ,1)
  90     tmp = (row-(r+1)).^2 + (col-(r+1)).^2;
  91     circ(row,col) = (tmp <= r^2);
  92   endfor
  93   circ = bwmorph(circ, 'remove');
  94 endfunction