--- a/inst/normxcorr2.m
+++ b/inst/normxcorr2.m
@@ -14,71 +14,127 @@
 ## this program; if not, see <http://www.gnu.org/licenses/>.
 
 ## -*- texinfo -*-
-## @deftypefn {Function File} {} normxcorr2 (@var{a}, @var{b})
-## Compute the 2D cross-correlation coefficient of matrices @var{a} and @var{b}.
+## @deftypefn {Function File} {} normxcorr2 (@var{template}, @var{img})
+## Compute normalized cross-correlation.
+##
+## Returns the the cross-correlation coefficient of matrices @var{template}
+## and @var{img}, a matrix of the same size as @var{img} with values ranging
+## between -1 and 1.
+##
+## Normalized correlation is mostly used for template matching, finding an
+## object or pattern, @var{template}, withing an image @var{img}.  Higher
+## values on the output show their locations, even in the presence of noise.
 ##
+## @group
+## @example
+## img = randi (255, 600, 400);
+## template = imnoise (img(100:150, 300:320), "gaussian");
+## cc = normxcorr2 (template, img);
+## [r, c] = find (cc == max (cc(:)))
+## @result{r} 150
+## @result{c} 320
+## @end example
+## @end group
 ##
-## @seealso{xcorr2, conv2, corr2, xcorr}
+## Despite the function name, this function will accept input with an arbitrary
+## number of dimensions.
+##
+## @seealso{conv2, convn, corr2, xcorr, xcorr2}
 ## @end deftypefn
 
+## Author: Benjamin Eltzner <b.eltzner@gmx.de>
+
 function c = normxcorr2 (a, b)
   if (nargin != 2)
-    print_usage;
-  endif
-  if (ndims (a) != 2 || ndims (b) != 2)
-    error ("normxcorr2: input matrices must have only 2 dimensions");
+    print_usage ();
   endif
 
-  ## compute cross correlation coefficient
-  [ma,na] = size(a);
-  [mb,nb] = size(b);
-  if (ma > mb || na > nb)
-    warning ("Template is larger than image.\nArguments may be accidentally interchanged.");
+  ## If this happens, it is probably a mistake
+  if (ndims (a) > ndims (b) || any (postpad (size (a), ndims (b)) > size (b)))
+    warning ("normxcorr2: TEMPLATE larger than IMG. Arguments may be swapped.");
   endif
 
-  a = double (a);
-  b = double (b);
-  a = a .- mean2(a);
-  b = b .- mean2(b);
-  c = conv2 (b, conj (a (ma:-1:1, na:-1:1)));
-  b = conv2 (b.^2, ones (size (a))) .- conv2 (b, ones (size (a))).^2 ./ (ma*na);
+  a = double (a) - mean (a(:));
+  b = double (b) - mean (b(:));
+
+  a1 = ones (size (a));
+  ar = reshape (a(end:-1:1), size (a));
+
+  c = convn (b, conj (ar));
+  b = convn (b.^2, a1) .- convn (b, a1).^2 ./ (prod (size (a)));
   a = sumsq (a(:));
-  c(:,:) = c(:,:) ./ sqrt (b(:,:) * a);
-  c(isnan(c)) = 0;
+  c = reshape (c ./ sqrt (b * a), size (c));
+
+  c(isnan (c)) = 0;
 endfunction
 
-%!test # basic usage
-%!shared a, b, c, row_shift, col_shift, a_dev1, b_dev1, a_dev2, b_dev2
+%!function offsets = get_max_offsets (c)
+%!  l = find (c == max (c(:)));
+%!  offsets = nthargout (1:ndims (c), @ind2sub, size (c), l);
+%!endfunction
+
+## test basic usage
+%!test
 %! row_shift = 18;
 %! col_shift = 20;
 %! a = randi (255, 30, 30);
 %! b = a(row_shift-10:row_shift, col_shift-7:col_shift);
 %! c = normxcorr2 (b, a);
-%!assert (nthargout ([1 2], @find, c == max (c(:))), {row_shift, col_shift}); # should return exact coordinates
-%! m = rand (size (b)) > 0.5;
-%! b(m) = b(m) * 0.95;
-%! b(!m) = b(!m) * 1.05;
+%! ## should return exact coordinates
+%! assert (get_max_offsets (c), {row_shift col_shift});
+%!
+%! ## Even with some small noise, should return exact coordinates
+%! b = imnoise (b, "gaussian");
 %! c = normxcorr2 (b, a);
-%!assert (nthargout ([1 2], @find, c == max (c(:))), {row_shift, col_shift}); # even with some small noise, should return exact coordinates
-%!test # coeff of autocorrelation must be same as negative of correlation by additive inverse
+%! assert (get_max_offsets (c), {row_shift col_shift});
+
+## The value for a "perfect" match should be 1. However, machine precision
+## creeps in most of the times.
+%!test
+%! a = rand (10, 10);
+%! c = normxcorr2 (a(5:7, 6:9), a);
+%! assert (c(7, 9), 1, eps*2);
+
+## coeff of autocorrelation must be same as negative of correlation
+## by additive inverse
+%!test
 %! a = 10 * randn (100, 100);
 %! auto = normxcorr2 (a, a);
 %! add_in = normxcorr2 (a, -a);
 %! assert (auto, -add_in);
-%!test # normalized correlation should be independent of scaling and shifting up to rounding errors
+
+## Normalized correlation should be independent of scaling and shifting
+## up to rounding errors
+%!test
 %! a = 10 * randn (50, 50);
 %! b = 10 * randn (100, 100);
-%! scale = 0;
-%! while (scale == 0)
-%! scale = 100 * rand();
-%! endwhile
-%! assert (max (max (normxcorr2 (scale*a,b) .- normxcorr2 (a,b))) < 1e-10);
-%! assert (max (max (normxcorr2 (a,scale*b) .- normxcorr2 (a,b))) < 1e-10);
+%! do
+%!   scale = 100 * rand ();
+%! until (scale != 0)
+%!
+%! assert (max ((normxcorr2 (scale*a,b) .- normxcorr2 (a,b))(:)), 0, 1e-10);
+%! assert (max ((normxcorr2 (a,scale*b) .- normxcorr2 (a,b))(:)), 0, 1e-10);
+%!
 %! a_shift1 = a .+ scale * ones (size (a));
 %! b_shift1 = b .+ scale * ones (size (b));
 %! a_shift2 = a .- scale * ones (size (a));
 %! b_shift2 = b .- scale * ones (size (b));
-%! assert (max (max (normxcorr2 (a_shift1,b) .- normxcorr2 (a,b))) < 1e-10);
-%! assert (max (max (normxcorr2 (a,b_shift1) .- normxcorr2 (a,b))) < 1e-10);
-%! assert (max (max (normxcorr2 (a_shift2,b) .- normxcorr2 (a,b))) < 1e-10);
-%! assert (max (max (normxcorr2 (a,b_shift2) .- normxcorr2 (a,b))) < 1e-10);
+%! assert (max ((normxcorr2 (a_shift1,b) .- normxcorr2 (a,b))(:)), 0, 1e-10);
+%! assert (max ((normxcorr2 (a,b_shift1) .- normxcorr2 (a,b))(:)), 0, 1e-10);
+%! assert (max ((normxcorr2 (a_shift2,b) .- normxcorr2 (a,b))(:)), 0, 1e-10);
+%! assert (max ((normxcorr2 (a,b_shift2) .- normxcorr2 (a,b))(:)), 0, 1e-10);
+
+## test n dimensional input
+%!test
+%! a = randi (100, 15, 15, 15);
+%! c = normxcorr2 (a(5:10, 2:6, 3:7), a);
+%! assert (get_max_offsets (c), {10 6 7});
+%!
+%! a = randi (100, 15, 15, 15);
+%! c = normxcorr2 (a(5:10, 2:6, 1:1), a);
+%! assert (get_max_offsets (c), {10 6 1});
+
+%!warning <swapped> normxcorr2 (rand (20), rand (5));
+%!error normxcorr2 (rand (5));
+%!error normxcorr2 (rand (5), rand (20), 2);
+