1## Copyright (C) 2013 Carnë Draug <carandraug@octave.org> 2## 3## This program is free software; you can redistribute it and/or modify it under 4## the terms of the GNU General Public License as published by the Free Software 5## Foundation; either version 3 of the License, or (at your option) any later 6## version. 7## 8## This program is distributed in the hope that it will be useful, but WITHOUT 9## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 10## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 11## details. 12## 13## You should have received a copy of the GNU General Public License along with 14## this program; if not, see <http://www.gnu.org/licenses/>. 15 16## -*- texinfo -*- 17## @deftypefn {Function File} {} nlfilter (@var{A}, @var{block_size}, @var{func}) 18## @deftypefnx {Function File} {} nlfilter (@var{A}, @var{block_size}, @var{func}, @dots{}) 19## @deftypefnx {Function File} {} nlfilter (@var{A}, "indexed", @dots{}) 20## Process matrix in sliding blocks with user-supplied function. 21## 22## Executes the function @var{fun} on each sliding block of 23## size @var{block_size}, taken from the matrix @var{A}. Both the matrix 24## @var{A}, and the block can have any number of dimensions. This function 25## is specially useful to perform sliding/moving window functions such as 26## moving average. 27## 28## The output will have the same dimensions @var{A}, each one of its values 29## corresponding to the processing of a block centered at the same 30## coordinates in @var{A}, with @var{A} being padded with zeros for the 31## borders (see below for indexed images). In case any side of the block 32## is of even length, the center is considered at indices 33## @code{floor ([@var{block_size}/2] + 1)}. 34## 35## The argument @var{func} must be a function handle that takes matrices of 36## size @var{block_size} as input and returns a single scalar. Any extra 37## input arguments to @code{nlfilter} are passed to @var{func} after the 38## block matrix. 39## 40## If @var{A} is an indexed image, the second argument should be the 41## string @qcode{"indexed"} so that any required padding is done correctly 42## as done by @code{im2col}. 43## 44## @emph{Note}: if @var{func} is a column compression function, i.e., it acts 45## along a column to return a single value, consider using @code{colfilt} 46## which usually performs faster. If @var{func} makes use of the colon 47## operator to select all elements in the block, e.g., if @var{func} looks 48## anything like @code{@@(x) sum (x(:))}, it is a good indication that 49## @code{colfilt} should be used. In addition, many sliding block operations 50## have their own specific implementations (see help text of @code{colfilt} 51## for a list). 52## 53## @seealso{blockproc, col2im, colfilt, im2col} 54## @end deftypefn 55 56function B = nlfilter (A, varargin) 57 58 ## Input check 59 [p, block_size, padval] = im2col_check ("nlfilter", nargin, A, varargin{:}); 60 61 if (nargin < p) 62 print_usage (); 63 endif 64 func = varargin{p++}; 65 if (! isa (func, "function_handle")) 66 error ("nlfilter: FUNC must be a function handle"); 67 elseif (! isscalar (func (ones (block_size, class (A)), varargin{p:nargin-1}))) 68 error ("nlfilter: FUNC must take a matrix of size BLOCK_SIZE and return a scalar"); 69 endif 70 ## Remaining params are parameters to fun. We need to place them into 71 ## a separate variable, we can't varargin{p:nargin-1} inside the anonymous 72 ## function because nargin will have a different value there 73 extra_args = varargin(p:nargin -1); 74 75 ## Pad image as required 76 padded = pad_for_sliding_filter (A, block_size, padval); 77 78 ## Get the blocks (one per column) 79 blks = im2col (padded, block_size, "sliding"); 80 81 ## New function that reshapes the array into a block before 82 ## passing it to the actual user supplied function 83 blk_func = @(x, z) func (reshape (blks(x:z), block_size), extra_args{:}); 84 85 ## Perform the filtering 86 blk_step = 1:(rows (blks)):(rows (blks) * columns (blks)); 87 B = arrayfun (blk_func, blk_step, blk_step + rows (blks) -1); 88 89 ## Back into its original shape 90 B = reshape (B, size (A)); 91 92endfunction 93 94%!demo 95%! ## creates a "wide" diagonal (although it can be performed more 96%! ## efficiently with "imdilate (A, true (3))") 97%! nlfilter (eye (10), [3 3], @(x) any (x(:) > 0)) 98 99%!assert (nlfilter (eye (4), [2 3], @(x) sum (x(:))), 100%! [2 2 1 0 101%! 1 2 2 1 102%! 0 1 2 2 103%! 0 0 1 1]); 104%!assert (nlfilter (eye (4), "indexed", [2 3], @(x) sum (x(:))), 105%! [4 2 1 2 106%! 3 2 2 3 107%! 2 1 2 4 108%! 4 3 4 5]); 109%!assert (nlfilter (eye (4), "indexed", [2 3], @(x, y) sum (x(:)) == y, 2), 110%! logical ([0 1 0 1 111%! 0 1 1 0 112%! 1 0 1 0 113%! 0 0 0 0])); 114 115## Check uint8 and uint16 padding (only the padding since the class of 116## the output is dependent on the function and sum() always returns double) 117%!assert (nlfilter (uint8 (eye (4)), "indexed", [2 3], @(x) sum (x(:))), 118%! [2 2 1 0 119%! 1 2 2 1 120%! 0 1 2 2 121%! 0 0 1 1]); 122%!assert (nlfilter (int16 (eye (4)), "indexed", [2 3], @(x) sum (x(:))), 123%! [4 2 1 2 124%! 3 2 2 3 125%! 2 1 2 4 126%! 4 3 4 5]); 127 128## Check if function class is preserved 129%!assert (nlfilter (uint8 (eye (4)), "indexed", [2 3], @(x) int8 (sum (x(:)))), 130%! int8 ([2 2 1 0 131%! 1 2 2 1 132%! 0 1 2 2 133%! 0 0 1 1])); 134 135%!test 136%! ## Effect of out of border elements. 137%! expected = [ 138%! 0.5 6.0 6.0 0.5 0 139%! 5.5 10.5 13.5 10.5 4.0 140%! 6.5 12.5 13.5 13.5 1.5 141%! 10.5 12.5 15.5 11.0 1.0 142%! 5.0 10.5 6.0 1.0 0 143%! ]; 144%! assert (nlfilter (magic (5), [3 4], @(x) median (x(:))), expected) 145 146 147%!test 148%! ## The center pixel of a sliding window when its length is even 149%! ## sized is ceil ((size (NHOOD) +1) /2) 150%! expected = [ 151%! 24 24 24 16 16 152%! 24 24 24 22 22 153%! 23 23 22 22 22 154%! 25 25 25 25 22 155%! 25 25 25 25 21 156%! ]; 157%! assert (nlfilter (magic (5), [3 4], @(x) max (x(:))), expected) 158 159## nlfilter and imdilate use a different center pixel when the 160## window/nhood have an even sized length. So to have imdilate behave 161## like nlfilter, we need to flip those dimensions. 162%!function dilated = imdilate_like_nlfilter (im, nhood) 163%! even_nhood_dims = find (mod (size (nhood), 2) == 0); 164%! for i = 1:even_nhood_dims 165%! im = flip (im, i); 166%! endfor 167%! dilated = imdilate (im, nhood); 168%! for i = 1:even_nhood_dims 169%! dilated = flip (dilated, i); 170%! endfor 171%!endfunction 172 173## Test N-dimensional 174%!test 175%! a = randi (65535, 20, 20, 20, "uint16"); 176%! ## extra dimensions on matrix only 177%! assert (nlfilter (a, [5 5], @(x) max(x(:))), imdilate (a, ones (5))) 178%! ## extra dimensions on both matrix and block 179%! assert (nlfilter (a, [5 5 5], @(x) max(x(:))), imdilate (a, ones ([5 5 5]))) 180%! ## extra dimensions and padding 181%! assert (nlfilter (a, [3 7], @(x) max(x(:))), imdilate (a, ones ([3 7]))) 182%! assert (nlfilter (a, [3 7 3], @(x) max(x(:))), imdilate (a, ones ([3 7 3]))) 183 184%!test 185%! a = randi (65535, 15, 15, 4, 8, 3, "uint16"); 186%! assert (nlfilter (a, [3 4 7 5], @(x) max(x(:))), 187%! imdilate_like_nlfilter (a, ones ([3 4 7 5]))) 188 189## Just to make sure it's not a bug in imdilate 190%!test 191%! a = randi (65535, 15, 15, 4, 3, 8, "uint16"); 192%! ord = ordfiltn (a, 3, ones ([3 7 3 1 5])); 193%! assert (nlfilter (a, [3 7 3 1 5], @(x) sort (x(:))(3)), ord) 194%! assert (nlfilter (a, [3 7 3 1 5], @(x, y) sort (x(:))(y), 3), ord) 195