octave_packages/image-1.0.15/nlfilter.m

   1 ## Copyright (C) 2004 Josep Mones i Teixidor
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 2 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {@var{B} = } nlfilter (@var{A}, [@var{m},@var{n}], @var{fun})
  18 ## @deftypefnx {Function File} {@var{B} = } nlfilter (@var{A}, [@var{m},@var{n}], @var{fun}, ...)
  19 ## @deftypefnx {Function File} {@var{B} = } nlfilter (@var{A},'indexed', ...)
  20 ## Processes image in sliding blocks using user-supplied function.
  21 ##
  22 ## @code{B=nlfilter(A,[m,n],fun)} passes sliding @var{m}-by-@var{n}
  23 ## blocks to user-supplied function @var{fun}. A block is build for
  24 ## every pixel in @var{A}, such as it is centered within the block.
  25 ## @var{fun} must return a scalar, and it is used to create matrix
  26 ## @var{B}. @var{nlfilter} pads the @var{m}-by-@var{n} block at the
  27 ## edges if necessary.
  28 ##
  29 ## Center of block is taken at ceil([@var{m},@var{n}]/2).
  30 ##
  31 ## @code{B=nlfilter(A,[m,n],fun,...)} behaves as described above but
  32 ## passes extra parameters to function @var{fun}.
  33 ##
  34 ## @code{B=nlfilter(A,'indexed',...)} assumes that @var{A} is an indexed
  35 ## image, so it pads the image using proper value: 0 for uint8 and
  36 ## uint16 images and 1 for double images. Keep in mind that if 'indexed'
  37 ## is not specified padding is always done using 0.
  38 ##
  39 ## @seealso{colfilt,blkproc,inline}
  40 ## @end deftypefn
  41
  42 ## Author:  Josep Mones i Teixidor <jmones@puntbarra.com>
  43
  44 function B = nlfilter(A, varargin)
  45   if(nargin<3)
  46     error("nlfilter: invalid number of parameters.");
  47   endif
  48
  49   ## check 'indexed' presence
  50   indexed=false;
  51   p=1;
  52   if(ischar(varargin{1}) && strcmp(varargin{1}, "indexed"))
  53     indexed=true;
  54     p+=1;
  55     if(isa(A,"uint8") || isa(A,"uint16"))
  56         padval=0;
  57     else
  58       padval=1;
  59     endif
  60   else
  61     padval=0;
  62   endif
  63
  64   ## check [m,n]
  65   if(!isvector(varargin{p}))
  66     error("nlfilter: expected [m,n] but param is not a vector.");
  67   endif
  68   if(length(varargin{p})!=2)
  69     error("nlfilter: expected [m,n] but param has wrong length.");
  70   endif
  71   sblk=varargin{p}(:);
  72   p+=1;
  73
  74   ## check fun
  75   ## TODO: add proper checks for this one
  76   if(nargin<p)
  77     error("nlfilter: required parameters haven't been supplied.");
  78   endif
  79   fun=varargin{p};
  80
  81   ## remaining params are params to fun
  82   ## extra params are p+1:nargin-1
  83
  84   ## We take an easy approach... feel free to optimize it (coding this
  85   ## in C++ would be a great idea).
  86
  87   ## Calculate center of block
  88   c=ceil(sblk/2);
  89
  90   ## Pre-padding
  91   prepad=c-ones(2,1);
  92
  93   ## Post-padding
  94   postpad=sblk-c;
  95
  96   ## Save A size
  97   as=size(A);
  98
  99   ## Pad data
 100   if(all(prepad==postpad))
 101     if(any(prepad>0))
 102       A=padarray(A,prepad,padval,'both');
 103     endif
 104   else
 105     if(any(prepad>0))
 106       A=padarray(A,prepad,padval,'pre');
 107     endif
 108     if(any(postpad>0))
 109       A=padarray(A,postpad,padval,'post');
 110     endif
 111   endif
 112
 113   ## calc end offsets
 114   me=postpad(1)+prepad(1);
 115   ne=postpad(2)+prepad(2);
 116
 117   ## We concatenate everything to preserve fun return type
 118   for i=1:as(1)
 119     r=feval(fun,A(i:i+me,1:1+ne),varargin{p+1:nargin-1});
 120     for j=2:as(2)
 121       r=horzcat(r,feval(fun,A(i:i+me,j:j+ne),varargin{p+1:nargin-1}));
 122     endfor
 123     if(i==1)
 124       B=r;
 125     else
 126       B=vertcat(B,r);
 127     endif
 128   endfor
 129
 130 endfunction
 131
 132 %!demo
 133 %! nlfilter(eye(10),[3,3],inline("any(x(:)>0)","x"))
 134 %! # creates a "wide" diagonal
 135
 136 %!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]);
 137 %!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]);
 138 %!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);
 139
 140 % Check uint8 and uint16 padding
 141 %!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]);
 142 %!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]);
 143
 144 % Check if function class is preserved
 145 %!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]));
 146
 147
 148
 149 %
 150 % $Log$
 151 % Revision 1.4  2007/03/23 16:14:37  adb014
 152 % Update the FSF address
 153 %
 154 % Revision 1.3  2007/01/04 23:50:47  hauberg
 155 % Put seealso before end deftypefn
 156 %
 157 % Revision 1.2  2007/01/04 23:37:54  hauberg
 158 % Minor changes in help text
 159 %
 160 % Revision 1.1  2006/08/20 12:59:35  hauberg
 161 % Changed the structure to match the package system
 162 %
 163 % Revision 1.5  2005/09/08 02:00:17  pkienzle
 164 % [for Bill Denney] isstr -> ischar
 165 %
 166 % Revision 1.4  2004/11/15 16:04:20  pkienzle
 167 % Fix tests for functions which return boolean matrices
 168 %
 169 % Revision 1.3  2004/09/03 13:28:32  jmones
 170 % Corrected behaviour for int* and uint* types
 171 %
 172 % Revision 1.2  2004/08/15 19:43:11  jmones
 173 % corrected a typo in doc
 174 %
 175 % Revision 1.1  2004/08/15 19:42:14  jmones
 176 % nlfilter: Processes image in siliding blocks using user-supplied function
 177 %
 178 %