# HG changeset patch # User Carnë Draug # Date 1395097714 0 # Node ID 0ed283edee13a5c562b51521a457d5d996fcf022 # Parent e9c18bff13be86a0d067969c2a3dcfc405edb0b2 normxcorr2: expand for N dimensions, new tests and documentation. * normxcorr2.m: use convn plus other changes so it accepts matrices or images with an arbitrary number of dimensions. Expand documentation with examples, and add tests for N dimensions, input checking, and precision. diff --git a/inst/normxcorr2.m b/inst/normxcorr2.m --- a/inst/normxcorr2.m +++ b/inst/normxcorr2.m @@ -14,71 +14,127 @@ ## this program; if not, see . ## -*- 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 + 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 normxcorr2 (rand (20), rand (5)); +%!error normxcorr2 (rand (5)); +%!error normxcorr2 (rand (5), rand (20), 2); +