--- a/COPYING
+++ b/COPYING
@@ -43,7 +43,7 @@
 inst/bwperim.m                          GPLv3+
 inst/bwselect.m                         GPLv3+
 inst/col2im.m                           GPLv3+
-inst/colfilt.m                          public domain
+inst/colfilt.m                          GPLv3+
 inst/colorgradient.m                    public domain
 inst/conndef.m                          GPLv3+
 inst/corr2.m                            GPLv3+

--- a/NEWS
+++ b/NEWS
@@ -126,6 +126,7 @@
 
       bestblk
       col2im
+      colfilt
       im2col
       nlfilter
 

--- a/inst/colfilt.m
+++ b/inst/colfilt.m
@@ -1,120 +1,181 @@
-## Author: Paul Kienzle <pkienzle@users.sf.net>
-## This program is granted to the public domain.
+## 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/>.
 
 ## -*- texinfo -*-
-## @deftypefn  {Function File} {} colfilt (@var{A}, [@var{r} @var{c}], [@var{m} @var{n}], "sliding", @var{f}, @dots{})
-## @deftypefnx {Function File} {} colfilt (@var{A}, [@var{r} @var{c}], [@var{m} @var{n}], "sliding", @var{f}, @dots{})
-## @deftypefnx {Function File} {} colfilt (@var{A}, [@var{r} @var{c}], [@var{m} @var{n}], "sliding", @var{f}, @dots{})
-## Apply filter to matrix blocks
+## @deftypefn  {Function File} {} colfilt (@var{A}, @var{block_size}, @var{block_type}, @var{func})
+## @deftypefnx {Function File} {} colfilt (@var{A}, @var{block_size}, @var{subsize}, @var{block_type}, @var{func}, @dots{})
+## @deftypefnx {Function File} {} colfilt (@var{A}, "indexed", @dots{})
+## @deftypefnx {Function File} {} colfilt (@dots{}, @var{func}, @var{extra_args}, @dots{})
+## Apply function to matrix blocks
+##
+## Executes the function @var{func} on blocks 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.
+##
+## The different blocks are organized into a matrix, each block as a
+## single column, and passed as the first to the function handle
+## @var{func}.  Any input arguments to @code{colfilt} after @var{func}
+## are passed to @var{func} after the blocks matrix.
+##
+## Blocks can be of two different types as defined by the string @var{block_type}:
 ##
-##   For each @var{r} x @var{c} overlapping subblock of @var{A}, add a column in matrix @var{C}
-##   @var{f}(@var{C},@dots{}) should return a row vector which is then reshaped into a
-##   a matrix of size @var{A} and returned. @var{A} is processed in chunks of size @var{m} x @var{n}.
+## @table @qcode
+## @item "distinct"
+## Each block is completely distinct from the other, with no overlapping
+## elements.  @var{func} must return a matrix of exactly the same size as
+## its input.
+##
+## @item "sliding"
+## Each possible block of size @var{block_size} inside @var{A} is used.
+## @var{func} should act along the column dimension (be a column
+## compression function) and return a vector of length equal to the
+## number of columns of its input.
+##
+## @end table
+##
+## The optional argument @var{subsize} divides @var{A} into smaller pieces
+## before generating the matrices with one block per column in order to
+## save memory.  It is currently only accepted for @sc{Matlab} compatibility.
 ##
-## The present version requires that [@var{m}, @var{n}] divide size(@var{A}), but for
-## compatibility it should work even if [@var{m}, @var{n}] does not divide @var{A}. Use
-## the following instead:
-## @example
-## [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);
-## @end example
+## If @var{A} is an indexed image, the second argument should be the
+## string @qcode{"indexed"} so that any required padding is done correctly.
+## The padding value will be 0 except for indexed images of class uint8
+## and uint16.
+##
+## This function is mostly useful to apply moving or sliding window filter
+## when @var{block_type} is "sliding".  However, for many cases, specialized
+## functions perform much faster.  For the following common cases, consider
+## the suggested alternatives;
+##
+## @table @asis
+## @item moving average
+## A moving average filter is equivalent to convolve with a matrix
+## of @code{1/@var{N}} sized @var{block_size}, where @var{N} is the total
+## number of elements in a block.  Use
+## @code{convn (@var{A}, (1/@var{N}) * ones (@var{block_size}) *, "same")}
+##
+## @item maximum or minimum
+## This is the equivalent to a dilation and erosion.  Use @code{imdilate} or
+## @code{imerode}.
 ##
-## @seealso{bestblk, blockproc, col2im, im2col}
+## @item any or all
+## Same as dilation and erosion but with logical input.  Use @code{imdilate}
+## or @code{imerode} with @code{logical (@var{A})}.
+##
+## @item median
+## Use @code{medfilt2} if @var{A} is only 2 dimensional, and @code{ordfiltn}
+## with the @code{floor (prod (@var{N}/ 2)} th element, where @var{N} is the
+## total number of elements in a block (add 1 if it is an even number).
+##
+## @item sort
+## Use @code{ordfiltn}.
+##
+## @item standard deviation
+## Use @code{stdfilt}.
+##
+## @item sum
+## Use a matrix of 1 to perform convolution,
+## @code{convn (@var{A}, ones (@var{block_size}), "same")}
+##
+## @end table
+##
+## @seealso{bestblk, blockproc, col2im, im2col, nlfilter}
 ## @end deftypefn
 
-## FIXME: distinct is not yet implemented
-## @deftypefn {Function File} {} colfilt (@var{A}, [@var{r}, @var{c}], [@var{m}, @var{n}], 'distinct', @var{f},@dots{})
-##   For each @var{r} x @var{c} non-overlapping subblock of @var{A}, add a column in matrix @var{C}
-##   @var{f}(@var{C},@dots{}) should return a matrix of size @var{C} each column of which is
-##   placed back into the subblock from whence it came. @var{A} is processed
-##   in chunks of size @var{m} x @var{n}.
+function B = colfilt (A, varargin)
+
+  ## Input check
+  if (nargin < 4)
+    print_usage ();
+  endif
+  [p, block_size, padval] = im2col_check ("colfilt", nargin, A, varargin{:});
 
-function B = colfilt(A,filtsize,blksize,blktype,f,varargin)
-   ## Input checking
-   real_nargin = nargin - length(varargin);
-   if (real_nargin < 4)
-     error("colfilt: not enough input arguments");
-   endif
-   if (ischar(blksize))
-     varargin = {f, varargin{:}};
-     f = blktype;
-     blktype = blksize;
-     blksize = size(A);
-   elseif (real_nargin < 5)
-     error("colfilt: not enough input arguments");
-   endif
-   if (!ismatrix(A) || ndims(A) != 2)
-     error("colfilt: first input argument must be a matrix");
-   endif
-   if (!isvector(filtsize) || numel(filtsize) != 2)
-     error("colfilt: second input argument must be a 2-vector");
-   endif
- 
-   ## Compute!
-   [m,n] = size(A);
-   r = filtsize(1);
-   c = filtsize(2);
-   mblock = blksize(1);
-   nblock = blksize(2);
+  subsize = size (A);
+  if (numel (varargin) < p)
+    print_usage ();
+  elseif (isnumeric (varargin{p}) && isvector (varargin{p}))
+    subsize = varargin{p++};
+    subsize = postpad (subsize, ndims (A), 1);
+    subsize = min (subsize, size (A));
+    subsize = max (subsize, block_size);
+  endif
 
-   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(:);
+  ## We need at least 2 more arguments (block type and function)
+  if (numel (varargin) < p +1)
+    print_usage ();
+  endif
+
+  ## Next one needs to be block type
+  block_type = varargin{p++};
+  if (! ischar (block_type))
+    error ("colfilt: BLOCK_TYPE must be a string");
+  endif
+
+  ## followed by the function
+  func = varargin{p++};
+  if (! isa (func, "function_handle"))
+    error ("colfilt: FUNC must be a function handle");
+  endif
 
-     # throw away old A to save memory.
-     B=A; clear A;  
+  ## anything after this are extra arguments to func
+  extra_args = varargin(p:end);
 
-     # 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;
+  switch (tolower (block_type))
+    case "sliding"
+      ## Function must return a single vector, one element per column,
+      ## i.e., should act along the elements of each column.
 
-     # 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,varargin{:});
-         B(1+i*mblock:(i+1)*mblock, 1+j*nblock:(j+1)*nblock) = tmp;
-       end
-     end
+      ## TODO for some large blocks, we may easily try to create matrix
+      ##      too large for Octave (see im2col documentation about the
+      ##      size). May be a good idea to split it into smaller images
+      ##      even if subsize is that large, so that we never go above
+      ##      sizemax ().
+      ##      However, this can be tricky. After splitting the image in
+      ##      smaller blocks, they can't be distinct, some parts need
+      ##      to overlap otherwise when we put them back together, we'll
+      ##      introduce many artifacts.
 
-   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,varargin{:});
-   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,varargin{:});
-     B = reshape(B,size(A));
-   endswitch
+      padded = pad_for_sliding_filter (A, block_size, padval);
+      cols = im2col (padded, block_size, "sliding");
+      B = col2im (func (cols, extra_args{:}), block_size, size (padded), "sliding");
+
+    case "distinct"
+      ## Function must return a matrix with the same number of elements
+      ## as its input, specially the same number of rows.
+
+      ## One of the options of this function is to do the block
+      ## processing already from big blocks from the original matrix
+      ## in order to save memory. While this may make sense with
+      ## sliding blocks, not so much here since cols will have the same
+      ## size as A, and so will B.
+      cols = im2col (A, block_size, "distinct");
+      B = col2im (func (cols, extra_args{:}), block_size, size (A), "distinct");
+
+    otherwise
+      error ("colfilt: invalid BLOCK_TYPE `%s'.", block_type);
+  endswitch
+
 endfunction
 
+%!demo
+%! ## Perform moving average filter with a 4x4 window
+%! A = magic (12)
+%! colfilt (A, [4 4], "sliding", @mean)
+
 %!test
-%! A = reshape(1:36,6,6);
-%! assert(colfilt(A,[2,2],[3,3],'sliding','sum'), conv2(A,ones(2),'same'));
+%! A = reshape (1:36, [6 6]);
+%! assert (colfilt (A, [2 2], [3 3], "sliding", @sum),
+%!         conv2 (A, ones (2), "same"));
+

--- a/inst/nlfilter.m
+++ b/inst/nlfilter.m
@@ -41,6 +41,15 @@
 ## string @qcode{"indexed"} so that any required padding is done correctly
 ## as done by @code{im2col}.
 ##
+## @emph{Note}: if @var{func} is a column compression function, i.e., it acts
+## along a column to return a single value, consider using @code{colfilt}
+## which usually performs faster.  If @var{func} makes use of the colon
+## operator to select all elements in the block, e.g., if @var{func} looks
+## anything like @code{@(x) sum (x(:))}, it is a good indication that
+## @code{colfilt} should be used.  In addition, many sliding block operations
+## have their own specific implementations (see help text of @code{colfilt}
+## for a list).
+##
 ## @seealso{blockproc, col2im, colfilt, im2col}
 ## @end deftypefn