Mercurial > hg > octave-image
annotate inst/immaximas.m @ 821:72cd72a3bff6
bestblk: complete rewrite to suppport multi-dimensional matrices.
* bestblk.m: implement support for any number of dimensions. Vectorize
code (loop only over dimensions now). Add tests.
* NEWS: new section for functions now supporting N-dimensional matrices.
| author | Carnë Draug <carandraug@octave.org> |
|---|---|
| date | Fri, 01 Nov 2013 05:05:28 +0000 |
| parents | a0ec22eef32c |
| children |
| rev | line source |
|---|---|
| 292 | 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/ | |
|
662
a0ec22eef32c
immaximas: fix license by using proper MIT license text after private conversation with copyright owner Peter Kovesi
carandraug
parents:
651
diff
changeset
|
5 ## |
| 292 | 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 | |
|
662
a0ec22eef32c
immaximas: fix license by using proper MIT license text after private conversation with copyright owner Peter Kovesi
carandraug
parents:
651
diff
changeset
|
8 ## in the Software without restriction, including without limitation the rights |
|
a0ec22eef32c
immaximas: fix license by using proper MIT license text after private conversation with copyright owner Peter Kovesi
carandraug
parents:
651
diff
changeset
|
9 ## to use, copy, modify, merge, publish, distribute, sublicense, and/or sell |
|
a0ec22eef32c
immaximas: fix license by using proper MIT license text after private conversation with copyright owner Peter Kovesi
carandraug
parents:
651
diff
changeset
|
10 ## copies of the Software, and to permit persons to whom the Software is |
|
a0ec22eef32c
immaximas: fix license by using proper MIT license text after private conversation with copyright owner Peter Kovesi
carandraug
parents:
651
diff
changeset
|
11 ## furnished to do so, subject to the following conditions: |
|
a0ec22eef32c
immaximas: fix license by using proper MIT license text after private conversation with copyright owner Peter Kovesi
carandraug
parents:
651
diff
changeset
|
12 ## |
|
a0ec22eef32c
immaximas: fix license by using proper MIT license text after private conversation with copyright owner Peter Kovesi
carandraug
parents:
651
diff
changeset
|
13 ## The above copyright notice and this permission notice shall be included in |
|
a0ec22eef32c
immaximas: fix license by using proper MIT license text after private conversation with copyright owner Peter Kovesi
carandraug
parents:
651
diff
changeset
|
14 ## all copies or substantial portions of the Software. |
|
a0ec22eef32c
immaximas: fix license by using proper MIT license text after private conversation with copyright owner Peter Kovesi
carandraug
parents:
651
diff
changeset
|
15 ## |
|
a0ec22eef32c
immaximas: fix license by using proper MIT license text after private conversation with copyright owner Peter Kovesi
carandraug
parents:
651
diff
changeset
|
16 ## The software is provided "as is", without warranty of any kind, express or |
|
a0ec22eef32c
immaximas: fix license by using proper MIT license text after private conversation with copyright owner Peter Kovesi
carandraug
parents:
651
diff
changeset
|
17 ## implied, including but not limited to the warranties of merchantability, |
|
a0ec22eef32c
immaximas: fix license by using proper MIT license text after private conversation with copyright owner Peter Kovesi
carandraug
parents:
651
diff
changeset
|
18 ## fitness for a particular purpose and noninfringement. In no event shall the |
|
a0ec22eef32c
immaximas: fix license by using proper MIT license text after private conversation with copyright owner Peter Kovesi
carandraug
parents:
651
diff
changeset
|
19 ## authors or copyright holders be liable for any claim, damages or other |
|
a0ec22eef32c
immaximas: fix license by using proper MIT license text after private conversation with copyright owner Peter Kovesi
carandraug
parents:
651
diff
changeset
|
20 ## liability, whether in an action of contract, tort or otherwise, arising from, |
|
a0ec22eef32c
immaximas: fix license by using proper MIT license text after private conversation with copyright owner Peter Kovesi
carandraug
parents:
651
diff
changeset
|
21 ## out of or in connection with the software or the use or other dealings in the |
|
a0ec22eef32c
immaximas: fix license by using proper MIT license text after private conversation with copyright owner Peter Kovesi
carandraug
parents:
651
diff
changeset
|
22 ## software. |
| 292 | 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}) | |
| 651 | 32 ## @deftypefnx{Function File} {[@var{r}, @var{c}, @dots{}] =} immaximas (@dots{}) |
| 33 ## @deftypefnx{Function File} {[@dots{}, @var{val}] =} immaximas (@dots{}) | |
| 292 | 34 ## Finds local spatial maximas of the given image. A local spatial maxima is |
| 35 ## defined as an image point with a value that is larger than all neighbouring | |
| 36 ## values in a square region of width 2*@var{radius}+1. By default @var{radius} | |
| 37 ## is 1, such that a 3 by 3 neighbourhood is searched. If the @var{thresh} input | |
| 38 ## argument is supplied, only local maximas with a value greater than @var{thresh} | |
| 39 ## are retained. | |
| 40 ## | |
| 41 ## The output vectors @var{r} and @var{c} contain the row-column coordinates | |
| 42 ## of the local maximas. The actual values are computed to sub-pixel precision | |
| 302 | 43 ## by fitting a parabola to the data around the pixel. If @var{im} is |
| 44 ## @math{N}-dimensional, then @math{N} vectors will be returned. | |
| 45 ## | |
| 46 ## If @var{im} is @math{N}-dimensional, and @math{N}+1 outputs are requested, | |
| 47 ## then the last output will contain the image values at the maximas. Currently | |
| 48 ## this value is not interpolated. | |
| 49 ## | |
| 50 ## @seealso{ordfilt2, ordfiltn} | |
| 292 | 51 ## @end deftypefn |
| 52 | |
| 302 | 53 function varargout = immaximas(im, radius, thresh) |
| 292 | 54 ## Check input |
| 55 if (nargin == 0) | |
| 56 error("immaximas: not enough input arguments"); | |
| 57 endif | |
| 58 if (nargin <= 1 || isempty(radius)) | |
| 59 radius = 1; | |
| 60 endif | |
| 61 if (nargin <= 2) | |
| 62 thresh = []; | |
| 63 endif | |
| 302 | 64 if (!ismatrix(im)) |
| 65 error("immaximas: first input argument must be an array"); | |
| 292 | 66 endif |
| 302 | 67 if (!isscalar(radius)) |
| 292 | 68 error("immaximas: second input argument must be a scalar or an empty matrix"); |
| 69 endif | |
| 70 if (!isscalar(thresh) && !isempty(thresh)) | |
| 71 error("immaximas: third input argument must be a scalar or an empty matrix"); | |
| 72 endif | |
| 73 | |
| 74 ## Find local maximas | |
| 302 | 75 nd = ndims(im); |
| 292 | 76 s = size(im); |
| 77 sze = 2*radius+1; | |
| 302 | 78 mx = ordfiltn(im, sze^nd, ones(repmat(sze,1, nd), "logical"), "reflect"); |
| 79 mx2 = ordfiltn(im, sze^nd-1, ones(repmat(sze,1, nd), "logical"), "reflect"); | |
| 292 | 80 |
| 302 | 81 # Find maxima, threshold |
| 82 immx = (im == mx) & (im != mx2); | |
| 292 | 83 if (!isempty(thresh)) |
| 84 immx &= (im>thresh); | |
| 85 endif | |
| 86 | |
| 87 ## Find local maximas and fit parabolas locally | |
| 302 | 88 ind = find(immx); |
| 89 [sub{1:nd}] = ind2sub(s, ind); | |
| 90 if (!isempty(ind)) | |
| 292 | 91 w = 1; # Width that we look out on each side of the feature point to fit a local parabola |
| 355 | 92 ws = w*cumprod([1; s(:)]); |
| 93 | |
| 302 | 94 ## We fit a parabola to the points in each dimension |
| 95 for d = 1:nd | |
| 96 ## Indices of points above, below, left and right of feature point | |
| 355 | 97 indminus1 = max(ind-ws(d), 1); |
| 98 indplus1 = min(ind+ws(d), numel(immx)); | |
| 292 | 99 |
| 302 | 100 ## Solve quadratic |
| 101 c = im(ind); | |
| 102 a = (im(indminus1) + im(indplus1))/2 - c; | |
| 103 b = a + c - im(indminus1); | |
| 104 shift = -w*b./(2*a); # Maxima of quadradic | |
| 105 | |
| 106 ## Move point | |
| 107 sub{d} += shift; | |
| 108 endfor | |
| 109 endif | |
| 110 | |
| 111 ## Output | |
| 112 varargout(1:nd) = sub(1:nd); | |
| 113 if (nargout > nd) | |
| 114 varargout{nd+1} = im(ind); | |
| 292 | 115 endif |
| 116 endfunction |