Mercurial > hg > octave-image
changeset 21:ea70f05f08b7
image block processing function
author | pkienzle |
---|---|
date | Fri, 10 May 2002 18:12:34 +0000 |
parents | 4de573bbaf4b |
children | e749d76cd428 |
files | colfilt.m |
diffstat | 1 files changed, 92 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/colfilt.m @@ -0,0 +1,92 @@ +## COLFILT Apply filter to matrix blocks +## colfilt(A,[r c],[m n],'sliding',f,...) +## For each r x c overlapping subblock of A, add a column in matrix C +## f(C,...) should return a row vector which is then reshaped into a +## a matrix of size A and returned. A is processed in chunks of size m x n. +## colfilt(A,[r c],[m n],'distinct',f,...) +## For each r x c non-overlapping subblock of A, add a column in matrix C +## f(C,...) should return a matrix of size C each column of which is +## placed back into the subblock from whence it came. A is processed +## in chunks of size m x n. +## +## The present version requires [m n], but for compatibility it should +## be optional. Use colfilt(A,[r c],size(A),...) +## +## The present version requires that [m n] divide size(A), but for +## compatibility it should work even if [m n] does not divide A. Use +## the following instead: +## [r c] = size(A); +## padA = zeros (m*ceil(r/m),n*ceil(c/n)); +## padA(1:r,1:c) = A; +## B = colfilt(padA,...); +## B = B(1:r,1:c); +## +## The present version does not handle 'distinct' + +## This software is granted to the public domain +function B = colfilt(A,filtsize,blksize,blktype,f,...) + + [m,n]=size(A); + r = filtsize(1); + c = filtsize(2); + mblock = blksize(1); + nblock = blksize(2); + + switch blktype + case 'sliding' + # pad with zeros + padm = (m+r-1); + padn = (n+c-1); + padA = zeros(padm, padn); + padA([1:m]+floor((r-1)/2),[1:n]+floor((c-1)/2)) = A; + padA = padA(:); + + # throw away old A to save memory. + B=A; clear A; + + # build the index vector + colidx = [0:r-1]'*ones(1,c) + padm*ones(r,1)*[0:c-1]; + offset = [1:mblock]'*ones(1,nblock) + padm*ones(mblock,1)*[0:nblock-1]; + idx = colidx(:)*ones(1,mblock*nblock) + ones(r*c,1)*offset(:)'; + clear colidx offset; + + # process the matrix, one block at a time + idxA = zeros(r*c,mblock*nblock); + tmp = zeros(mblock,nblock); + for i = 0:m/mblock-1 + for j = 0:n/nblock-1 + idxA(:) = padA(idx + (i*mblock + padm*j*nblock)); + tmp(:) = feval(f,idxA,all_va_args); + B(1+i*mblock:(i+1)*mblock, 1+j*nblock:(j+1)*nblock) = tmp; + end + end + + case 'old-sliding' # processes the whole matrix at a time + padA = zeros(m+r-1,n+c-1); + padA([1:m]+floor(r/2),[1:n]+floor(c/2)) = A; + [padm,padn] = size(padA); + colidx = [0:r-1]'*ones(1,c) + padm*ones(r,1)*[0:c-1]; + offset = [1:m]'*ones(1,n) + padm*ones(m,1)*[0:n-1]; + idx = colidx(:)*ones(1,m*n) + ones(r*c,1)*offset(:)'; + idxA = zeros(r*c,m*n); + idxA(:) = padA(:)(idx); + B = zeros(size(A)); + B(:) = feval(f,idxA,all_va_args); + case 'old-distinct' # processes the whole matrix at a time + if (r*floor(m/r) != m || c*floor(n/c) != n) + error("colfilt expected blocks to exactly fill A"); + endif + colidx = [0:r-1]'*ones(1,c) + m*ones(r,1)*[0:c-1]; + offset = [1:r:m]'*ones(1,n/c) + m*ones(m/r,1)*[0:c:n-1]; + idx =colidx(:)*ones(1,m*n/r/c) + ones(r*c,1)*offset(:)'; + idxA = zeros(r*c,m*n/r/c); + idxA(:) = A(:)(idx); + B = zeros(prod(size(A)),1); + B(idx) = feval(f,idxA,all_va_args); + B = reshape(B,size(A)); + endswitch +endfunction + +%!test +%! A = reshape(1:36,6,6); +%! assert(colfilt(A,[2,2],[3,3],'sliding','sum'), conv2(A,ones(2),'same');