2 ## @deftypefn {Function File} colfilt(@var{A}, [@var{r}, @var{c}], [@var{m}, @var{n}], 'sliding', @var{f},...)
3 ## Apply filter to matrix blocks
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}.
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:
19 ## padA = zeros (m*ceil(r/m),n*ceil(c/n));
21 ## B = colfilt(padA,...);
25 ## The present version does not handle 'distinct'
28 ## This software is granted to the public domain
29 ## Author: Paul Kienzle <pkienzle@users.sf.net>
31 function B = colfilt(A,filtsize,blksize,blktype,f,varargin)
33 real_nargin = nargin - length(varargin);
35 error("colfilt: not enough input arguments");
38 varargin = {f, varargin{:}};
42 elseif (real_nargin < 5)
43 error("colfilt: not enough input arguments");
45 if (!ismatrix(A) || ndims(A) != 2)
46 error("colfilt: first input argument must be a matrix");
48 if (!isvector(filtsize) || numel(filtsize) != 2)
49 error("colfilt: second input argument must be a 2-vector");
64 padA = zeros(padm, padn);
65 padA([1:m]+floor((r-1)/2),[1:n]+floor((c-1)/2)) = A;
68 # throw away old A to save memory.
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(:)';
77 # process the matrix, one block at a time
78 idxA = zeros(r*c,mblock*nblock);
79 tmp = zeros(mblock,nblock);
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;
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);
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");
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);
108 B = zeros(prod(size(A)),1);
109 B(idx) = feval(f,idxA,varargin{:});
110 B = reshape(B,size(A));
115 %! A = reshape(1:36,6,6);
116 %! assert(colfilt(A,[2,2],[3,3],'sliding','sum'), conv2(A,ones(2),'same'));