Rewrite normalization algorithm.

This is a major revision that removes the `method` parameter of `match_template` and uses a new normalization method.

Note that the example result is different with this new normalization.

doc/examples/plot_template.py

@@ -29,8 +29,7 @@ for x, y in target_positions:
     image[x:x+size, y:y+size] = target
 image += randn(400, 400)*2
 
-# Match the template.
-result = match_template(image, target, method='norm-corr')
+result = match_template(image, target)
 
 found_positions = peak_local_max(result)
 

skimage/feature/_template.pyx

@@ -1,4 +1,34 @@
-"""template.py - Template matching
+"""
+Template matching using normalized cross-correlation.
+
+We use fast normalized cross-correlation algorithm (see [1]_ and [2]_) to
+compute match probability. This algorithm calculates the normalized
+cross-correlation of an image, `I`, with a template `T` according to the
+following equation::
+
+                sum{ I(x, y) [T(x, y) - <T>] }
+    -------------------------------------------------------
+    sqrt(sum{ [I(x, y) - <I>]^2 } sum{ [T(x, y) - <T>]^2 })
+
+where `<T>` is the average of the template, and `<I>` is the average of the
+image *coincident with the template*, and sums are over the template and the
+image window coincident with the template. Note that the numerator is simply
+the cross-correlation of the image and the zero-mean template.
+
+To speed up calculations, we use summed-area tables (a.k.a. integral images) to
+quickly calculate sums of image windows inside the loop. This step relies on
+the following relation (see Eq. 10 of [1])::
+
+    sum{ [I(x, y) - <I>]^2 } =
+        sum{ I^2(x, y) } - [sum{ I(x, y) }]^2 / N_x N_y
+
+(Without this relation, you would need to subtract each image-window mean from
+the image window *before* squaring.)
+
+.. [1] Briechle and Hanebeck, "Template Matching using Fast Normalized
+       Cross Correlation", Proceedings of the SPIE (2001).
+.. [2] J. P. Lewis, "Fast Normalized Cross-Correlation", Industrial Light and
+       Magic.
 """
 import cython
 cimport numpy as np
@@ -54,36 +84,25 @@ cdef float sum_integral(np.ndarray[float, ndim=2,  mode="c"] sat,
 
 @cython.boundscheck(False)
 def match_template(np.ndarray[float, ndim=2, mode="c"] image,
-                   np.ndarray[float, ndim=2, mode="c"] template,
-                   str method):
-    cdef np.ndarray[float, ndim=2] result
-    cdef np.ndarray[float, ndim=2,  mode="c"] integral_sum
-    cdef np.ndarray[float, ndim=2,  mode="c"] integral_sqr
+                   np.ndarray[float, ndim=2, mode="c"] template):
+    cdef np.ndarray[float, ndim=2, mode="c"] result
+    cdef np.ndarray[float, ndim=2, mode="c"] integral_sum
+    cdef np.ndarray[float, ndim=2, mode="c"] integral_sqr
     cdef float template_mean = np.mean(template)
-    cdef float template_norm
+    cdef float template_ssd
     cdef float inv_area
 
+    integral_sum = integral.integral_image(image)
+    integral_sqr = integral.integral_image(image**2)
+
+    template -= template_mean
+    template_ssd = np.sum(template**2)
+    # use inversed area for accuracy
+    inv_area = 1.0 / (template.shape[0] * template.shape[1])
+
     # when `dtype=float` is used, ascontiguousarray returns ``double``.
     result = np.ascontiguousarray(fftconvolve(image, np.fliplr(template),
                                               mode="valid"), dtype=np.float32)
-    if method == 'norm-coeff':
-        integral_sum = integral.integral_image(image)
-    integral_sqr = integral.integral_image(image**2)
-
-    # use inversed area for accuracy
-    inv_area = 1.0 / (template.shape[0] * template.shape[1])
-
-    if method == 'norm-corr':
-        # calculate template norm according to the following:
-        # variance = 1/K Sum[(x_k - mean) ** 2]
-        #          = 1/K Sum[x_k ** 2] - mean ** 2
-        #template_norm = sqrt((np.std(template) ** 2 +
-                              #template_mean ** 2)) / sqrt(inv_area)
-        # TODO: check equation for template_norm.
-        # The above normalization factor is equivalent to the second-moment.
-        template_norm = sqrt(np.sum(template**2))
-    else:
-        template_norm = sqrt((template_mean ** 2)) / sqrt(inv_area)
 
     cdef int i, j
     cdef float num, den, window_sqr_sum, window_mean_sqr, window_sum,
@@ -96,15 +115,11 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image,
             i_end = i + template.shape[0] - 1
             j_end = j + template.shape[1] - 1
 
-            window_mean_sqr = 0
-            if method == 'norm-coeff':
-                window_sum = sum_integral(integral_sum, i, j, i_end, j_end)
-                window_mean_sqr = window_sum * window_sum * inv_area
-                num -= window_sum * template_mean
-
+            window_sum = sum_integral(integral_sum, i, j, i_end, j_end)
+            window_mean_sqr = window_sum * window_sum * inv_area
             window_sqr_sum = sum_integral(integral_sqr, i, j, i_end, j_end)
+            den = sqrt((window_sqr_sum - window_mean_sqr) * template_ssd)
 
-            den = sqrt(window_sqr_sum - window_mean_sqr) * template_norm
             # enforce some limits
             if fabs(num) < den:
                 num /= den

skimage/feature/template.py

@@ -6,11 +6,12 @@ import _template
 from skimage.util.dtype import _convert
 
 
-def match_template(image, template, method='norm-coeff', pad_output=True):
-    """Finds a template in an image using normalized correlation.
+def match_template(image, template, pad_output=True):
+    """Match a template to an image using normalized correlation.
 
-    TODO: The output is currently smaller than the input image due to
-          cropping at the boundaries equal to the template width.
+    The output is an array with values between -1.0 and 1.0, which correspond
+    to the probability that the template's *origin* (i.e. its top-left
+    corner) is found at that position.
 
     Parameters
     ----------
@@ -18,23 +19,6 @@ def match_template(image, template, method='norm-coeff', pad_output=True):
         Image to process.
     template : array_like
         Template to locate.
-    method : str
-        The correlation method used in scanning.
-        T represents the template, I the image and R the result.
-        All sums are done over X = 0..w-1 and Y = 0..h-1 of the template.
-        'norm-coeff':
-            R(x, y) = Sum[T(X, Y) * I(x + X, y + Y)] / N
-            N = sqrt(Sum[T(X, Y)**2] * Sum[I(x + X, y + Y)**2])
-        'norm-corr':
-            R(x,y) = Sum[T'(X, Y) * I'(x + X, y + Y)] / N
-            N = sqrt(Sum[T'(X, Y)**2] * Sum[I'(x + X, y + Y)**2])
-
-            where:
-
-            T'(x, y) = T(X, Y) - mean(T)
-            I'(x + X, y + Y) = I(x + X, y + Y) - mean[I(X', Y')]
-            mean[I(X', Y')] = mean of image region under the template.
-
     pad_output : bool
         If True, pad output array to be the same size as the input image.
         Otherwise, the output is an array with shape `(M - m + 1, N - n + 1)`
@@ -43,18 +27,15 @@ def match_template(image, template, method='norm-coeff', pad_output=True):
     Returns
     -------
     output : ndarray
-        Correlation results between 0.0 and 1.0, which correspond to the match
-        probability when the template's *origin* (i.e. its top-left corner) is
-        placed at that position. The bottom and right edges of `output` are
-        truncated (`pad_output = False`) or zero-padded (`pad_output = True`),
-        since otherwise the template would extend beyond the image edges.
+        Correlation results between -1.0 and 1.0. The `output` is truncated
+        (`pad_output = False`) or zero-padded (`pad_output = True`) at the
+        bottom and right edges, where the template would otherwise extend
+        beyond the image edges.
 
     """
-    if method not in ('norm-corr', 'norm-coeff'):
-        raise ValueError("Unknown template method: %s" % method)
     image = _convert(image, np.float32)
     template = _convert(template, np.float32)
-    result = _template.match_template(image, template, method)
+    result = _template.match_template(image, template)
 
     if pad_output:
         h, w = result.shape

skimage/feature/tests/test_template.py

@@ -1,5 +1,4 @@
 import numpy as np
-from numpy.random import randn
 from numpy.testing import assert_array_almost_equal as assert_close
 
 from skimage.feature import match_template, peak_local_max
@@ -14,25 +13,25 @@ def test_template():
     target_positions = [(50, 50), (200, 200)]
     for x, y in target_positions:
         image[x:x + size, y:y + size] = target
-    image += randn(400, 400) * 2
+    np.random.seed(1)
+    image += np.random.randn(400, 400) * 2
 
-    for method in ["norm-corr", "norm-coeff"]:
-        result = match_template(image, target, method=method)
-        delta = 5
+    result = match_template(image, target)
+    delta = 5
 
-        positions = peak_local_max(result, min_distance=delta)
+    positions = peak_local_max(result, min_distance=delta)
 
-        if len(positions) > 2:
-            # Keep the two maximum peaks.
-            intensities = result[tuple(positions.T)]
-            i_maxsort = np.argsort(intensities)[::-1]
-            positions = positions[i_maxsort][:2]
+    if len(positions) > 2:
+        # Keep the two maximum peaks.
+        intensities = result[tuple(positions.T)]
+        i_maxsort = np.argsort(intensities)[::-1]
+        positions = positions[i_maxsort][:2]
 
-        # Sort so that order matches `target_positions`.
-        positions = positions[np.argsort(positions[:, 0])]
+    # Sort so that order matches `target_positions`.
+    positions = positions[np.argsort(positions[:, 0])]
 
-        for xy_target, xy in zip(target_positions, positions):
-            yield assert_close, xy, xy_target
+    for xy_target, xy in zip(target_positions, positions):
+        yield assert_close, xy, xy_target
 
 
 if __name__ == "__main__":