changeset 822:896f5b53db6b

Rewrite nlfilter for support of N-dimensional input. * nlfilter.m: complete rewrite to support N-dimensional images and blocks. Also make use of im2col rather than rewrite the block creation. * private/im2col_check.m: new function that performs that perform most of input checking for both im2col() and nlfilter(). * im2col.m: leave most input check to new private function im2col_check(), and add a few more tests for input checking. * private/pad_for_spatial_filter.m: rename to pad_for_sliding_filter for a better match of what it does (no longer only meant to prepare input for the __spatial_filter__ function). * private/pad_for_sliding_filter.m: previously named pad_for_spatial_filter(). Also changed one of its input argument to be size of block/domain, rather than the block/domain itself. * ordfiltn.m: fix call to the changed pad_for_spatial_filter function. * NEWS: add nlfilter to list of functions that now support N-dimensional input.
author Carnë Draug <carandraug@octave.org>
date Fri, 01 Nov 2013 19:24:40 +0000
parents 72cd72a3bff6
children 939332b5142e
files NEWS inst/im2col.m inst/nlfilter.m inst/ordfiltn.m inst/private/im2col_check.m inst/private/pad_for_sliding_filter.m inst/private/pad_for_spatial_filter.m
diffstat 6 files changed, 189 insertions(+), 156 deletions(-) [+]
line wrap: on
line diff
--- a/NEWS
+++ b/NEWS
@@ -127,6 +127,7 @@
       bestblk
       col2im
       im2col
+      nlfilter
 
  Summary of important user-visible changes for image 2.0.0 (2012/11/08):
 -------------------------------------------------------------------------
--- a/inst/im2col.m
+++ b/inst/im2col.m
@@ -69,54 +69,23 @@
 ## N-dimensional conversion of image blocks into colums.
 
 function B = im2col (A, varargin)
-  if (nargin < 2 || nargin > 4)
-    print_usage ();
-  elseif (! ismatrix (A) || ! (isnumeric (A) || islogical (A)))
-    error ("im2col: A must be a numeric of logical matrix");
-  endif
-  p = 1;  # varargin param being processsed
 
-  ## Defaults
-  block_type  = "sliding";
-  padval      = 0;
-  indexed     = false;
-
-  ## Check for 'indexed' presence
-  if (ischar (varargin{p}) && strcmpi (varargin{p}, "indexed"))
-    indexed = true;
-    if (nargin < 3)
-      print_usage ();
-    endif
-    ## We pad with value of 0 for indexed images of uint8 or uint16 class.
-    ## Indexed images of signed integer, or above uint16, are treated the
-    ## same as floating point (a value of 1 is index 1 on the colormap)
-    if (any (isa (A, {"uint8", "uint16"})))
-      padval = 0;
-    else
-      padval = 1;
-    endif
-    p++;
-  elseif (nargin > 3)
-    ## If we didn't have "indexed" but had 4 parameters there's an error
+  ## Input check
+  if (nargin > 4)
     print_usage ();
   endif
-
-  ## check [m,n]
-  block_size = varargin{p};
-  if (! isnumeric (block_size) || ! isvector (block_size) ||
-      any (block_size(:) < 1))
-    error ("im2col: BLOCK_SIZE must be a vector of positive elements.");
-  endif
-  block_size(end+1:ndims(A)) = 1; # expand singleton dimensions if required
-  block_size = block_size(:).';   # make sure it's a row vector
-  p++;
-
+  [p, block_size, padval] = im2col_check ("im2col", nargin, A, varargin{:});
   if (nargin > p)
     ## we have block_type param
     if (! ischar (varargin{p}))
-      error("im2col: invalid parameter block_type.");
+      error("im2col: BLOCK_TYPE must be a string");
     endif
-    block_type = varargin{p};
+    block_type = varargin{p++};
+  else
+    block_type = "sliding";
+  endif
+  if (nargin > p)
+    print_usage ();
   endif
 
   ## After all the input check, start the actual im2col. The idea is to
@@ -203,8 +172,11 @@
 %! a = rand (10);
 %! assert (im2col (a, [5 5]), im2col (a, "indexed", [5 5]))
 
+%!error <BLOCK_TYPE> im2col (rand (20), [2 5], 10)
 %!error <BLOCK_TYPE> im2col (rand (20), [2 5], "wrong_block_type")
 %!error <greater than> im2col (rand (10), [11 5], "sliding")
+%!error im2col (rand (10), [5 5], "sliding", 5)
+%!error im2col (rand (10), "indexed", [5 5], "sliding", 5)
 
 %!shared B, A, Bs, As, Ap, Bp0, Bp1, Bp0_3s
 %! v   = [1:10]';
--- a/inst/nlfilter.m
+++ b/inst/nlfilter.m
@@ -1,4 +1,4 @@
-## Copyright (C) 2004 Josep Mones i Teixidor <jmones@puntbarra.com>
+## Copyright (C) 2013 Carnë Draug <carandraug@octave.org>
 ##
 ## This program is free software; you can redistribute it and/or modify it under
 ## the terms of the GNU General Public License as published by the Free Software
@@ -14,130 +14,135 @@
 ## this program; if not, see <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {@var{B} = } nlfilter (@var{A}, [@var{m},@var{n}], @var{fun})
-## @deftypefnx {Function File} {@var{B} = } nlfilter (@var{A}, [@var{m},@var{n}], @var{fun}, @dots{})
-## @deftypefnx {Function File} {@var{B} = } nlfilter (@var{A},'indexed', @dots{})
-## Processes image in sliding blocks using user-supplied function.
+## @deftypefn  {Function File} {} nlfilter (@var{A}, @var{block_size}, @var{func})
+## @deftypefnx {Function File} {} nlfilter (@var{A}, @var{block_size}, @var{func}, @dots{})
+## @deftypefnx {Function File} {} nlfilter (@var{A}, "indexed", @dots{})
+## Process matrix in sliding blocks with user-supplied function.
+##
+## Executes the function @var{fun} on each sliding block of
+## size @var{block_size}, taken from the matrix @var{A}.  Both the matrix
+## @var{A}, and the block can have any number of dimensions.  This function
+## is specially useful to perform sliding/moving window functions such as
+## moving average.
 ##
-## @code{B=nlfilter(A,[m,n],fun)} passes sliding @var{m}-by-@var{n}
-## blocks to user-supplied function @var{fun}. A block is build for
-## every pixel in @var{A}, such as it is centered within the block.
-## @var{fun} must return a scalar, and it is used to create matrix
-## @var{B}. @var{nlfilter} pads the @var{m}-by-@var{n} block at the
-## edges if necessary.
-## 
-## Center of block is taken at ceil([@var{m},@var{n}]/2).
+## The output will have the same dimensions @var{A}, each one of its values
+## corresponding to the processing of a block centered at the same
+## coordinates in @var{A}, with @var{A} being padded with zeros for the
+## borders (see below for indexed images).  In case any side of the block
+## is of even length, the center is considered at indices
+## @code{floor ([@var{block_size}/2] + 1)}.
 ##
-## @code{B=nlfilter(A,[m,n],fun, @dots{})} behaves as described above but
-## passes extra parameters to function @var{fun}.
+## The argument @var{func} must be a function handle that takes matrices of
+## size @var{block_size} as input and returns a single scalar.  Any extra
+## input arguments to @code{nlfilter} are passed to @var{func} after the
+## block matrix.
 ##
-## @code{B=nlfilter(A,'indexed', @dots{})} assumes that @var{A} is an indexed
-## image, so it pads the image using proper value: 0 for uint8 and
-## uint16 images and 1 for double images. Keep in mind that if 'indexed'
-## is not specified padding is always done using 0.
+## If @var{A} is an indexed image, the second argument should be the
+## string @qcode{"indexed"} so that any required padding is done correctly
+## as done by @code{im2col}.
 ##
-## @seealso{colfilt,blkproc,inline}
+## @seealso{blockproc, col2im, colfilt, im2col}
 ## @end deftypefn
 
-function B = nlfilter(A, varargin)
-  if(nargin<3)
-    error("nlfilter: invalid number of parameters.");
-  endif
-  
-  ## check 'indexed' presence
-  indexed=false;
-  p=1;
-  if(ischar(varargin{1}) && strcmp(varargin{1}, "indexed"))
-    indexed=true;
-    p+=1;
-    if(isa(A,"uint8") || isa(A,"uint16"))
-      padval=0;
-    else
-      padval=1; 
-    endif
-  else
-    padval=0;
-  endif
+function B = nlfilter (A, varargin)
 
-  ## check [m,n]
-  if(!isvector(varargin{p}))
-    error("nlfilter: expected [m,n] but param is not a vector.");
-  endif
-  if(length(varargin{p})!=2)
-    error("nlfilter: expected [m,n] but param has wrong length.");
-  endif
-  sblk=varargin{p}(:);
-  p+=1;
-
-  ## check fun
-  ## TODO: add proper checks for this one
-  if(nargin<p)
-    error("nlfilter: required parameters haven't been supplied.");
-  endif
-  fun=varargin{p};
-  
-  ## remaining params are params to fun
-  ## extra params are p+1:nargin-1
+  ## Input check
+  [p, block_size, padval] = im2col_check ("nlfilter", nargin, A, varargin{:});
 
-  ## We take an easy approach... feel free to optimize it (coding this
-  ## in C++ would be a great idea).
-  
-  ## Calculate center of block
-  c=ceil(sblk/2);
-  
-  ## Pre-padding
-  prepad=c-ones(2,1);
-  
-  ## Post-padding
-  postpad=sblk-c;
-  
-  ## Save A size
-  as=size(A);
+  if (nargin < p)
+    print_usage ();
+  endif
+  func = varargin{p++};
+  if (! isa (func, "function_handle"))
+    error ("nlfilter: FUNC must be a function handle");
+  elseif (! isscalar (func (ones (block_size, class (A)), varargin{p:nargin-1})))
+    error ("nlfilter: FUNC must take a matrix of size BLOCK_SIZE and return a scalar");
+  endif
+  ## Remaining params are parameters to fun. We need to place them into
+  ## a separate variable, we can't varargin{p:nargin-1} inside the anonymous
+  ## function because nargin will have a different value there
+  extra_args = varargin(p:nargin -1);
 
-  ## Pad data
-  if(all(prepad==postpad))
-    if(any(prepad>0))
-      A=padarray(A,prepad,padval,'both');
-    endif
-  else
-    if(any(prepad>0))
-      A=padarray(A,prepad,padval,'pre');
-    endif
-    if(any(postpad>0))
-      A=padarray(A,postpad,padval,'post');
-    endif
-  endif
+  ## Pad image as required
+  padded = pad_for_sliding_filter (A, block_size, padval);
+
+  ## Get the blocks (one per column)
+  blks = im2col (padded, block_size, "sliding");
 
-  ## calc end offsets
-  me=postpad(1)+prepad(1);
-  ne=postpad(2)+prepad(2);
+  ## New function that reshapes the array into a block before
+  ## passing it to the actual user supplied function
+  blk_func = @(x, z) func (reshape (blks(x:z), block_size), extra_args{:});
 
-  ## We concatenate everything to preserve fun return type
-  for i=1:as(1)
-    r=feval(fun,A(i:i+me,1:1+ne),varargin{p+1:nargin-1});
-    for j=2:as(2)
-      r=horzcat(r,feval(fun,A(i:i+me,j:j+ne),varargin{p+1:nargin-1}));
-    endfor
-    if(i==1)
-      B=r;
-    else
-      B=vertcat(B,r);
-    endif
-  endfor
+  ## Perform the filtering
+  blk_step = 1:(rows (blks)):(rows (blks) * columns (blks));
+  B = arrayfun (blk_func, blk_step, blk_step + rows (blks) -1);
+
+  ## Back into its original shape
+  B = reshape (B, size (A));
 
 endfunction
 
 %!demo
-%! nlfilter(eye(10),[3,3],inline("any(x(:)>0)","x"))
-%! # creates a "wide" diagonal
+%! ## creates a "wide" diagonal (although it can be performed more
+%! ## efficiently with "imdilate (A, true (3))")
+%! nlfilter (eye (10), [3 3], @(x) any (x(:) > 0))
 
-%!assert(nlfilter(eye(4),[2,3],inline("sum(x(:))","x")),[2,2,1,0;1,2,2,1;0,1,2,2;0,0,1,1]);
-%!assert(nlfilter(eye(4),'indexed',[2,3],inline("sum(x(:))","x")),[4,2,1,2;3,2,2,3;2,1,2,4;4,3,4,5]);
-%!assert(nlfilter(eye(4),'indexed',[2,3],inline("sum(x(:))==y","x","y"),2),[0,1,0,1;0,1,1,0;1,0,1,0;0,0,0,0]!=0);
+%!assert (nlfilter (eye (4), [2 3], @(x) sum (x(:))),
+%!        [2 2 1 0
+%!         1 2 2 1
+%!         0 1 2 2
+%!         0 0 1 1]);
+%!assert (nlfilter (eye (4), "indexed", [2 3], @(x) sum (x(:))),
+%!        [4 2 1 2
+%!         3 2 2 3
+%!         2 1 2 4
+%!         4 3 4 5]);
+%!assert (nlfilter (eye (4), "indexed", [2 3], @(x, y) sum (x(:)) == y, 2),
+%!        logical ([0 1 0 1
+%!                  0 1 1 0
+%!                  1 0 1 0
+%!                  0 0 0 0]));
 
-% Check uint8 and uint16 padding
-%!assert(nlfilter(uint8(eye(4)),'indexed',[2,3],inline("sum(x(:))","x")),[2,2,1,0;1,2,2,1;0,1,2,2;0,0,1,1]);
-%!assert(nlfilter(uint16(eye(4)),'indexed',[2,3],inline("sum(x(:))","x")),[2,2,1,0;1,2,2,1;0,1,2,2;0,0,1,1]);
+## Check uint8 and uint16 padding (only the padding since the class of
+## the output is dependent on the function and sum() always returns double)
+%!assert (nlfilter (uint8 (eye (4)), "indexed", [2 3], @(x) sum (x(:))),
+%!        [2 2 1 0
+%!         1 2 2 1
+%!         0 1 2 2
+%!         0 0 1 1]);
+%!assert (nlfilter (int16 (eye (4)), "indexed", [2 3], @(x) sum (x(:))),
+%!        [4 2 1 2
+%!         3 2 2 3
+%!         2 1 2 4
+%!         4 3 4 5]);
+
+## Check if function class is preserved
+%!assert (nlfilter (uint8 (eye (4)), "indexed", [2 3], @(x) int8 (sum (x(:)))),
+%!        int8 ([2 2 1 0
+%!               1 2 2 1
+%!               0 1 2 2
+%!               0 0 1 1]));
 
-% Check if function class is preserved
-%!assert(nlfilter(uint8(eye(4)),'indexed',[2,3],inline("int8(sum(x(:)))","x")),int8([2,2,1,0;1,2,2,1;0,1,2,2;0,0,1,1]));
+## Test N-dimensional
+%!test
+%! a = randi (10, 20, 20, 20);
+%! ## extra dimensions on matrix only
+%! assert (nlfilter (a, [5 5], @(x) max(x(:))), imdilate (a, ones (5)))
+%! ## extra dimensions on both matrix and block
+%! assert (nlfilter (a, [5 5 5], @(x) max(x(:))), imdilate (a, ones ([5 5 5])))
+%! ## extra dimensions and padding
+%! assert (nlfilter (a, [3 7], @(x) max(x(:))), imdilate (a, ones ([3 7])))
+%! assert (nlfilter (a, [3 7 3], @(x) max(x(:))), imdilate (a, ones ([3 7 3])))
+
+## FIXME: this is failing due to a bug in imdilate
+%!test
+%! a = randi (10, 15, 15, 4, 8, 3);
+%! assert (nlfilter (a, [3 4 7 5], @(x) max(x(:))),
+%!         imdilate (a, ones ([3 4 7 5])))
+
+## Just to make sure it's not a bug in imdilate
+%!test
+%! a = randi (10, 15, 15, 4, 3, 8);
+%! ord = ordfiltn (a, 3, ones ([3 7 3 1 5]));
+%! assert (nlfilter (a, [3 7 3 1 5], @(x) sort (x(:))(3)), ord)
+%! assert (nlfilter (a, [3 7 3 1 5], @(x, y) sort (x(:))(y), 3), ord)
--- a/inst/ordfiltn.m
+++ b/inst/ordfiltn.m
@@ -77,7 +77,7 @@
     endif
   endfor
 
-  A = pad_for_spatial_filter (A, domain, padding);
+  A = pad_for_sliding_filter (A, size (domain), padding);
 
   ## Perform the filtering
   retval = __spatial_filtering__ (A, logical (domain), "ordered", S, nth);
new file mode 100644
--- /dev/null
+++ b/inst/private/im2col_check.m
@@ -0,0 +1,55 @@
+## Copyright (C) 2013 Carnë Draug <carandraug@octave.org>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+## A private function originally written for im2col but that can also
+## be used by colfilt and nlfilt, since the syntax is very similar.
+
+function [p, block_size, padval] = im2col_check (func, nargin, A, varargin)
+
+  if (nargin < 2)
+    print_usage (func);
+  elseif (! ismatrix (A) || ! (isnumeric (A) || islogical (A)))
+    error ("%s: A must be a numeric of logical matrix", func);
+  endif
+  p = 1;  # varargin param being processsed
+
+  padval = 0;
+
+  ## Check for 'indexed' presence
+  if (ischar (varargin{p}) && strcmpi (varargin{p}, "indexed"))
+    if (nargin < 3)
+      print_usage (func);
+    endif
+    ## We pad with value of 0 for indexed images of uint8 or uint16 class.
+    ## Indexed images of signed integer, or above uint16, are treated the
+    ## same as floating point (a value of 1 is index 1 on the colormap)
+    if (any (isa (A, {"uint8", "uint16"})))
+      padval = 0;
+    else
+      padval = 1;
+    endif
+    p++;
+  endif
+
+  ## check [m,n]
+  block_size = varargin{p++};
+  if (! isnumeric (block_size) || ! isvector (block_size) ||
+      any (block_size(:) < 1))
+    error ("%s: BLOCK_SIZE must be a vector of positive elements.", func);
+  endif
+  block_size(end+1:ndims(A)) = 1; # expand singleton dimensions if required
+  block_size = block_size(:).';   # make sure it's a row vector
+
+endfunction
rename from inst/private/pad_for_spatial_filter.m
rename to inst/private/pad_for_sliding_filter.m
--- a/inst/private/pad_for_spatial_filter.m
+++ b/inst/private/pad_for_sliding_filter.m
@@ -18,15 +18,15 @@
 ## __spatial_filtering__. It may be better to have __spatial_filtering__ do this
 ## part as well then
 
-function im = pad_for_spatial_filter (im, domain, padval)
+function im = pad_for_sliding_filter (im, window_size, padval)
 
-  if (ndims (im) != ndims (domain))
+  if (ndims (im) != numel (window_size))
     error ("image and domain must have the same number of dimensions")
   endif
 
-  pad  = floor (size (domain) / 2);
+  pad  = floor (window_size / 2);
   im   = padarray (im, pad, padval);
-  even = ! mod (size(domain), 2);
+  even = ! mod (window_size, 2);
 
   ## if one of the domain dimensions is even, its origin is must be
   ##     floor ([size(domain)/2] + 1)