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