octave_packages/image-1.0.15/immaximas.m

   1 ## Copyright (c) 2003-2005 Peter Kovesi
   2 ## School of Computer Science & Software Engineering
   3 ## The University of Western Australia
   4 ## http://www.csse.uwa.edu.au/
   5 ##
   6 ## Permission is hereby granted, free of charge, to any person obtaining a copy
   7 ## of this software and associated documentation files (the "Software"), to deal
   8 ## in the Software without restriction, subject to the following conditions:
   9 ##
  10 ## The above copyright notice and this permission notice shall be included in all
  11 ## copies or substantial portions of the Software.
  12 ##
  13 ## The Software is provided "as is", without warranty of any kind.
  14 ##
  15 ## I've made minor changes compared to the original 'nonmaxsuppts' function developed
  16 ## by Peter Kovesi. The original is available at
  17 ## http://www.csse.uwa.edu.au/~pk/research/matlabfns/Spatial/nonmaxsuppts.m
  18 ##    -- Søren Hauberg, 2008
  19
  20 ## -*- texinfo -*-
  21 ## @deftypefn {Function File} {[@var{r}, @var{c}] =} immaximas (@var{im}, @var{radius})
  22 ## @deftypefnx{Function File} {[@var{r}, @var{c}] =} immaximas (@var{im}, @var{radius}, @var{thresh})
  23 ## @deftypefnx{Function File} {[@var{r}, @var{c}, ...] =} immaximas (...)
  24 ## @deftypefnx{Function File} {[..., @var{val}] =} immaximas (...)
  25 ## Finds local spatial maximas of the given image. A local spatial maxima is
  26 ## defined as an image point with a value that is larger than all neighbouring
  27 ## values in a square region of width 2*@var{radius}+1. By default @var{radius}
  28 ## is 1, such that a 3 by 3 neighbourhood is searched. If the @var{thresh} input
  29 ## argument is supplied, only local maximas with a value greater than @var{thresh}
  30 ## are retained.
  31 ##
  32 ## The output vectors @var{r} and @var{c} contain the row-column coordinates
  33 ## of the local maximas. The actual values are computed to sub-pixel precision
  34 ## by fitting a parabola to the data around the pixel. If @var{im} is
  35 ## @math{N}-dimensional, then @math{N} vectors will be returned.
  36 ##
  37 ## If @var{im} is @math{N}-dimensional, and @math{N}+1 outputs are requested,
  38 ## then the last output will contain the image values at the maximas. Currently
  39 ## this value is not interpolated.
  40 ##
  41 ## @seealso{ordfilt2, ordfiltn}
  42 ## @end deftypefn
  43
  44 function varargout = immaximas(im, radius, thresh)
  45   ## Check input
  46   if (nargin == 0)
  47     error("immaximas: not enough input arguments");
  48   endif
  49   if (nargin <= 1 || isempty(radius))
  50     radius = 1;
  51   endif
  52   if (nargin <= 2)
  53     thresh = [];
  54   endif
  55   if (!ismatrix(im))
  56     error("immaximas: first input argument must be an array");
  57   endif
  58   if (!isscalar(radius))
  59     error("immaximas: second input argument must be a scalar or an empty matrix");
  60   endif
  61   if (!isscalar(thresh) && !isempty(thresh))
  62     error("immaximas: third input argument must be a scalar or an empty matrix");
  63   endif
  64
  65   ## Find local maximas
  66   nd = ndims(im);
  67   s = size(im);
  68   sze = 2*radius+1;
  69   mx  = ordfiltn(im, sze^nd,   ones(repmat(sze,1, nd), "logical"), "reflect");
  70   mx2 = ordfiltn(im, sze^nd-1, ones(repmat(sze,1, nd), "logical"), "reflect");
  71
  72   # Find maxima, threshold
  73   immx = (im == mx) & (im != mx2);
  74   if (!isempty(thresh))
  75     immx &= (im>thresh);
  76   endif
  77
  78   ## Find local maximas and fit parabolas locally
  79   ind = find(immx);
  80   [sub{1:nd}] = ind2sub(s, ind);
  81   if (!isempty(ind))
  82     w = 1; # Width that we look out on each side of the feature point to fit a local parabola
  83     ws = w*cumprod([1; s(:)]);
  84
  85     ## We fit a parabola to the points in each dimension
  86     for d = 1:nd
  87       ## Indices of points above, below, left and right of feature point
  88       indminus1 = max(ind-ws(d), 1);
  89       indplus1  = min(ind+ws(d), numel(immx));
  90
  91       ## Solve quadratic
  92       c = im(ind);
  93       a = (im(indminus1) + im(indplus1))/2 - c;
  94       b = a + c - im(indminus1);
  95       shift = -w*b./(2*a); # Maxima of quadradic
  96
  97       ## Move point
  98       sub{d} += shift;
  99     endfor
 100   endif
 101
 102   ## Output
 103   varargout(1:nd) = sub(1:nd);
 104   if (nargout > nd)
 105     varargout{nd+1} = im(ind);
 106   endif
 107 endfunction
 108