From 5682d27eb089108bb3c106e22ae6a7ce6baa0c91 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Thu, 1 Mar 2012 23:00:34 -0500 Subject: [PATCH] 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 | 3 +- skimage/feature/_template.pyx | 79 +++++++++++++++----------- skimage/feature/template.py | 39 ++++--------- skimage/feature/tests/test_template.py | 29 +++++----- 4 files changed, 72 insertions(+), 78 deletions(-) diff --git a/doc/examples/plot_template.py b/doc/examples/plot_template.py index fe6549f6..3aea93d9 100644 --- a/doc/examples/plot_template.py +++ b/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) diff --git a/skimage/feature/_template.pyx b/skimage/feature/_template.pyx index 9052fd3c..63ffa3d9 100644 --- a/skimage/feature/_template.pyx +++ b/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) - ] } + ------------------------------------------------------- + sqrt(sum{ [I(x, y) - ]^2 } sum{ [T(x, y) - ]^2 }) + +where `` is the average of the template, and `` 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) - ]^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 diff --git a/skimage/feature/template.py b/skimage/feature/template.py index 39713e43..3c85985f 100644 --- a/skimage/feature/template.py +++ b/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 diff --git a/skimage/feature/tests/test_template.py b/skimage/feature/tests/test_template.py index 0ec48f09..a3af9954 100644 --- a/skimage/feature/tests/test_template.py +++ b/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__":