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, including without limitation the rights
9## to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
10## copies of the Software, and to permit persons to whom the Software is
11## furnished to do so, subject to the following conditions:
12##
13## The above copyright notice and this permission notice shall be included in
14## all copies or substantial portions of the Software.
15##
16## The software is provided "as is", without warranty of any kind, express or
17## implied, including but not limited to the warranties of merchantability,
18## fitness for a particular purpose and noninfringement. In no event shall the
19## authors or copyright holders be liable for any claim, damages or other
20## liability, whether in an action of contract, tort or otherwise, arising from,
21## out of or in connection with the software or the use or other dealings in the
22## software.
23##
24## I've made minor changes compared to the original 'nonmaxsuppts' function developed
25## by Peter Kovesi. The original is available at
26## http://www.csse.uwa.edu.au/~pk/research/matlabfns/Spatial/nonmaxsuppts.m
27##    -- Søren Hauberg, 2008
28
29## -*- texinfo -*-
30## @deftypefn {Function File} {[@var{r}, @var{c}] =} immaximas (@var{im}, @var{radius})
31## @deftypefnx{Function File} {[@var{r}, @var{c}] =} immaximas (@var{im}, @var{radius}, @var{thresh})
32## @deftypefnx{Function File} {[@var{r}, @var{c}, @dots{}] =} immaximas (@dots{})
33## @deftypefnx{Function File} {[@dots{}, @var{val}] =} immaximas (@dots{})
34## Find local spatial maximas.
35##
36## Local spatial maximas should not be mistaken with regional maxima.
37## See @code{imregionalmax} for the later.
38##
39## A local spatial maxima is
40## defined as an image point with a value that is larger than all neighbouring
41## values in a square region of width 2*@var{radius}+1. By default @var{radius}
42## is 1, such that a 3 by 3 neighbourhood is searched. If the @var{thresh} input
43## argument is supplied, only local maximas with a value greater than @var{thresh}
44## are retained.
45##
46## The output vectors @var{r} and @var{c} contain the row-column coordinates
47## of the local maximas. The actual values are computed to sub-pixel precision
48## by fitting a parabola to the data around the pixel. If @var{im} is
49## @math{N}-dimensional, then @math{N} vectors will be returned.
50##
51## If @var{im} is @math{N}-dimensional, and @math{N}+1 outputs are requested,
52## then the last output will contain the image values at the maximas. Currently
53## this value is not interpolated.
54##
55## @seealso{imregionalmax, ordfilt2, ordfiltn}
56## @end deftypefn
57
58function varargout = immaximas(im, radius, thresh)
59  ## Check input
60  if (nargin == 0)
61    error("immaximas: not enough input arguments");
62  endif
63  if (nargin <= 1 || isempty(radius))
64    radius = 1;
65  endif
66  if (nargin <= 2)
67    thresh = [];
68  endif
69  if (! isnumeric (im))
70    error("immaximas: IM must be a numeric array");
71  endif
72  if (!isscalar(radius))
73    error("immaximas: second input argument must be a scalar or an empty matrix");
74  endif
75  if (!isscalar(thresh) && !isempty(thresh))
76    error("immaximas: third input argument must be a scalar or an empty matrix");
77  endif
78
79  ## Find local maximas
80  nd = ndims(im);
81  s = size(im);
82  sze = 2*radius+1;
83  mx  = ordfiltn(im, sze^nd,   ones(repmat(sze,1, nd), "logical"), "reflect");
84  mx2 = ordfiltn(im, sze^nd-1, ones(repmat(sze,1, nd), "logical"), "reflect");
85
86  # Find maxima, threshold
87  immx = (im == mx) & (im != mx2);
88  if (!isempty(thresh))
89    immx &= (im>thresh);
90  endif
91
92  ## Find local maximas and fit parabolas locally
93  ind = find(immx);
94  [sub{1:nd}] = ind2sub(s, ind);
95  if (!isempty(ind))
96    w = 1; # Width that we look out on each side of the feature point to fit a local parabola
97    ws = w*cumprod([1; s(:)]);
98
99    ## We fit a parabola to the points in each dimension
100    for d = 1:nd
101      ## Indices of points above, below, left and right of feature point
102      indminus1 = max(ind-ws(d), 1);
103      indplus1  = min(ind+ws(d), numel(immx));
104
105      ## Solve quadratic
106      c = im(ind);
107      a = (im(indminus1) + im(indplus1))/2 - c;
108      b = a + c - im(indminus1);
109      shift = -w*b./(2*a); # Maxima of quadratic
110
111      ## Move point
112      sub{d} += shift;
113    endfor
114  endif
115
116  ## Output
117  varargout(1:nd) = sub(1:nd);
118  if (nargout > nd)
119    varargout{nd+1} = im(ind);
120  endif
121endfunction
122