octave_packages/image-1.0.15/colfilt.m

   1 ## -*- texinfo -*-
   2 ## @deftypefn {Function File} colfilt(@var{A}, [@var{r}, @var{c}], [@var{m}, @var{n}], 'sliding', @var{f},...)
   3 ## Apply filter to matrix blocks
   4 ##
   5 ##   For each @var{r} x @var{c} overlapping subblock of @var{A}, add a column in matrix @var{C}
   6 ##   @var{f}(@var{C},...) should return a row vector which is then reshaped into a
   7 ##   a matrix of size @var{A} and returned. @var{A} is processed in chunks of size @var{m} x @var{n}.
   8 ## @deftypefnx{Function File} colfilt(@var{A}, [@var{r}, @var{c}], [@var{m}, @var{n}], 'distinct', @var{f},...)
   9 ##   For each @var{r} x @var{c} non-overlapping subblock of @var{A}, add a column in matrix @var{C}
  10 ##   @var{f}(@var{C},...) should return a matrix of size @var{C} each column of which is
  11 ##   placed back into the subblock from whence it came. @var{A} is processed
  12 ##   in chunks of size @var{m} x @var{n}.
  13 ##
  14 ## The present version requires that [@var{m}, @var{n}] divide size(@var{A}), but for
  15 ## compatibility it should work even if [@var{m}, @var{n}] does not divide @var{A}. Use
  16 ## the following instead:
  17 ## @example
  18 ## [r, c] = size(A);
  19 ## padA = zeros (m*ceil(r/m),n*ceil(c/n));
  20 ## padA(1:r,1:c) = A;
  21 ## B = colfilt(padA,...);
  22 ## B = B(1:r,1:c);
  23 ## @end example
  24 ##
  25 ## The present version does not handle 'distinct'
  26 ## @end deftypefn
  27
  28 ## This software is granted to the public domain
  29 ## Author: Paul Kienzle <pkienzle@users.sf.net>
  30
  31 function B = colfilt(A,filtsize,blksize,blktype,f,varargin)
  32    ## Input checking
  33    real_nargin = nargin - length(varargin);
  34    if (real_nargin < 4)
  35      error("colfilt: not enough input arguments");
  36    endif
  37    if (ischar(blksize))
  38      varargin = {f, varargin{:}};
  39      f = blktype;
  40      blktype = blksize;
  41      blksize = size(A);
  42    elseif (real_nargin < 5)
  43      error("colfilt: not enough input arguments");
  44    endif
  45    if (!ismatrix(A) || ndims(A) != 2)
  46      error("colfilt: first input argument must be a matrix");
  47    endif
  48    if (!isvector(filtsize) || numel(filtsize) != 2)
  49      error("colfilt: second input argument must be a 2-vector");
  50    endif
  51
  52    ## Compute!
  53    [m,n] = size(A);
  54    r = filtsize(1);
  55    c = filtsize(2);
  56    mblock = blksize(1);
  57    nblock = blksize(2);
  58
  59    switch blktype
  60    case 'sliding'
  61      # pad with zeros
  62      padm = (m+r-1);
  63      padn = (n+c-1);
  64      padA = zeros(padm, padn);
  65      padA([1:m]+floor((r-1)/2),[1:n]+floor((c-1)/2)) = A;
  66      padA = padA(:);
  67
  68      # throw away old A to save memory.
  69      B=A; clear A;
  70
  71      # build the index vector
  72      colidx = [0:r-1]'*ones(1,c) + padm*ones(r,1)*[0:c-1];
  73      offset = [1:mblock]'*ones(1,nblock) + padm*ones(mblock,1)*[0:nblock-1];
  74      idx = colidx(:)*ones(1,mblock*nblock) + ones(r*c,1)*offset(:)';
  75      clear colidx offset;
  76
  77      # process the matrix, one block at a time
  78      idxA = zeros(r*c,mblock*nblock);
  79      tmp = zeros(mblock,nblock);
  80      for i = 0:m/mblock-1
  81        for j = 0:n/nblock-1
  82          idxA(:) = padA(idx + (i*mblock + padm*j*nblock));
  83          tmp(:) = feval(f,idxA,varargin{:});
  84          B(1+i*mblock:(i+1)*mblock, 1+j*nblock:(j+1)*nblock) = tmp;
  85        end
  86      end
  87
  88    case 'old-sliding'  # processes the whole matrix at a time
  89      padA = zeros(m+r-1,n+c-1);
  90      padA([1:m]+floor(r/2),[1:n]+floor(c/2)) = A;
  91      [padm,padn] = size(padA);
  92      colidx = [0:r-1]'*ones(1,c) + padm*ones(r,1)*[0:c-1];
  93      offset = [1:m]'*ones(1,n) + padm*ones(m,1)*[0:n-1];
  94      idx = colidx(:)*ones(1,m*n) + ones(r*c,1)*offset(:)';
  95      idxA = zeros(r*c,m*n);
  96      idxA(:) = padA(:)(idx);
  97      B = zeros(size(A));
  98      B(:) = feval(f,idxA,varargin{:});
  99    case 'old-distinct' # processes the whole matrix at a time
 100      if (r*floor(m/r) != m || c*floor(n/c) != n)
 101         error("colfilt expected blocks to exactly fill A");
 102      endif
 103      colidx = [0:r-1]'*ones(1,c) + m*ones(r,1)*[0:c-1];
 104      offset = [1:r:m]'*ones(1,n/c) + m*ones(m/r,1)*[0:c:n-1];
 105      idx =colidx(:)*ones(1,m*n/r/c) + ones(r*c,1)*offset(:)';
 106      idxA = zeros(r*c,m*n/r/c);
 107      idxA(:) = A(:)(idx);
 108      B = zeros(prod(size(A)),1);
 109      B(idx) = feval(f,idxA,varargin{:});
 110      B = reshape(B,size(A));
 111    endswitch
 112 endfunction
 113
 114 %!test
 115 %! A = reshape(1:36,6,6);
 116 %! assert(colfilt(A,[2,2],[3,3],'sliding','sum'), conv2(A,ones(2),'same'));