--- a/NEWS
+++ b/NEWS
@@ -121,6 +121,13 @@
       bwperim
       padarray
 
+ ** The following functions now fully support matrices with an arbitrary
+    number of dimensions:
+
+      bestblk
+      col2im
+      im2col
+
  Summary of important user-visible changes for image 2.0.0 (2012/11/08):
 -------------------------------------------------------------------------
 

--- a/inst/bestblk.m
+++ b/inst/bestblk.m
@@ -1,4 +1,4 @@
-## Copyright (C) 2004 Josep Mones i Teixidor <jmones@puntbarra.com>
+## 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
@@ -14,72 +14,66 @@
 ## this program; if not, see <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {@var{siz} = } bestblk ([@var{m} @var{n}], @var{k})
-## @deftypefnx {Function File} {[@var{mb} @var{nb}] = } bestblk ([@var{m} @var{n}], @var{k})
-## Calculates the best size of block for block processing.
+## @deftypefn  {Function File} {@var{blk_siz} =} bestblk (@var{IMS})
+## @deftypefnx {Function File} {@var{blk_siz} =} bestblk (@var{IMS}, @var{max})
+## @deftypefnx {Function File} {[@var{Mb}, @var{Nb}, @dots{}] =} bestblk (@dots{})
+## Calculate block best size for block processing.
 ##
-## @code{siz=bestblk([m,n],k)} calculates the optimal block size for block
-## processing for a @var{m}-by-@var{n} image. @var{k} is the maximum
-## side dimension of the block. Its default value is 100. @var{siz} is a
-## row vector which contains row and column dimensions for the block.
+## Given a matrix of size @var{IMS}, calculates the largest size for distinct
+## blocks @var{blk_size}, that minimize padding and is smaller than or equal to
+## @var{k} (defaults to 100) 
 ##
-## @code{[mb,nb]=bestblk([m,n],k)} behaves as described above but
-## returns block dimensions to @var{mb} and @var{nb}.
+## The output @var{blk_siz] is a row vector for the block size.  If there are
+## multiple output arguments, the number of rows is assigned to the first
+## (@var{Mb}), and the number of columns to the second (@var{Nb}), etc.
 ##
-## @strong{Algorithm:}
+## To determine @var{block_size}, the following is performed for each
+## dimension:
 ##
-## For each dimension (@var{m} and @var{n}), it follows this algorithm:
-##
-## 1.- If dimension is less or equal than @var{k}, it returns the
+## @enumerate
+## @item
+## If dimension @var{IMS} is less or equal than @var{k}, it returns the
 ## dimension value.
 ##
-## 2.- If not then returns the value between
-## @code{round(min(dimension/10,k/2))} which minimizes padding.
+## @item
+## If not, find the highest value between @code{min (dimension/10, k/2)}
+## which minimizes padding.
 ##
+## @end enumerate
 ##
-## @seealso{blockproc}
+## @seealso{blockproc, col2im, im2col}
 ## @end deftypefn
 
 function [varargout] = bestblk (ims, k = 100)
-  if (nargin < 1 || nargin > 2 || nargout > 2)
-    print_usage;
-  elseif (!isvector (ims))
-    error("bestblk: first parameter is not a vector.");
-  endif
-  ims=ims(:);
-  if(length(ims)!=2)
-    error("bestblk: length of first parameter is not 2.");
+
+  if (nargin < 1 || nargin > 2)
+    print_usage ();
+  elseif (! isnumeric (ims) || ! isvector (ims) || any (ims(:) < 1))
+    error("bestblk: IMS must be a numeric vector of positive integers.");
+  elseif (numel (ims) < 2)
+    error ("bestblk: IMS must have at least 2 elements");
+  elseif (! isnumeric (k) || ! isscalar (k) || k < 1)
+    error ("bestblk: K must be a positive scalar");
   endif
 
-  mb=mi=ims(1);
-  p=mi;
-  if(mi>k)
-    for i=round(min(mi/10,k/2)):k
-      pt=rem(mi,i);
-      if(pt<p)
-        p=pt;
-        mb=i;
-      endif
-    endfor
-  endif
+  ims = floor (ims(:).');
+  k   = floor (k);
 
-  nb=ni=ims(2);
-  p=ni;
-  if(ni>k)
-    for i=round(min(ni/10,k/2)):k
-      pt=rem(ni,i);
-      if(pt<p)
-        p=pt;
-        nb=i;
-      endif
-    endfor
-  endif
+  out = zeros (size (ims));
+  for dim = 1:numel (ims)
+    if (ims(dim) <= k)
+      out(dim) = ims(dim);
+    else
+      possible = k:-1:min (ims(dim) /10, k /2);
+      [~, ind] = min (mod (-ims(dim), possible));
+      out(dim) = possible(ind);
+    endif
+  endfor
 
-  if(nargout<=1)
-    varargout{1}=[mb;nb];
+  if (nargout <= 1)
+    varargout{1} = out;
   else
-    varargout{1}=mb;
-    varargout{2}=nb;
+    varargout = mat2cell (out', ones (1, numel (out)));
   endif
 
 endfunction
@@ -88,5 +82,22 @@
 %! siz = bestblk ([200; 10], 50);
 %! disp (siz)
 
-%!assert(bestblk([300;100],150),[30;100]);
-%!assert(bestblk([256,128],17),[16;16]);
+%!error <numeric vector> bestblk ("string")
+%!error <positive scalar> bestblk ([100 200], "string")
+%!error <2 elements> bestblk ([100], 5)
+
+%!assert (bestblk ([ 10  12],   2), [  2   2]);
+%!assert (bestblk ([ 10  12],   3), [  2   3]);
+%!assert (bestblk ([300 100], 150), [150 100]);
+%!assert (bestblk ([256 128],  17), [ 16  16]);
+
+## make sure we really pick the highest one
+%!assert (bestblk ([ 17  17],   3), [  3   3]);
+
+## Test default
+%!assert (bestblk ([230 470]), bestblk ([230 470], 100))
+
+## Test N-dimensional
+%!assert (bestblk ([10 12 10], 3), [2 3 2]);
+%!assert (bestblk ([ 9 12  9], 3), [3 3 3]);
+%!assert (bestblk ([10 12 10 11], 5), [5 4 5 4]);