From 545bdab985a0c56f041b501601397e2586dfa694 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 4 Feb 2012 13:23:56 -0500 Subject: [PATCH] Refactor template matching. * Change Cython function to take names of correlation method instead of numbers representing the methods. * Use alternate formula for `template_norm` of 'norm-corr' method, but note that both formulas need to be checked for correctness. * Add note that `match_template` output has a different shape than the input image. This needs to be fixed before merging. * Change 'Sigma' to 'Sum' in docstring to avoid confusion with standard deviation. * Other minor changes for readability. --- skimage/feature/_template.pyx | 42 ++++++++++++++++++----------------- skimage/feature/template.py | 26 +++++++++++----------- 2 files changed, 35 insertions(+), 33 deletions(-) diff --git a/skimage/feature/_template.pyx b/skimage/feature/_template.pyx index d64edf04..44d78a47 100644 --- a/skimage/feature/_template.pyx +++ b/skimage/feature/_template.pyx @@ -55,30 +55,35 @@ 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, - int num_type): + str method): # convolve the image with template by frequency domain multiplication cdef np.ndarray[float, ndim=2] result # when `dtype=float` is used, ascontiguousarray returns ``double``. result = np.ascontiguousarray(fftconvolve(image, np.fliplr(template), mode="valid"), dtype=np.float32) + # calculate squared integral images used for normalization cdef np.ndarray[float, ndim=2, mode="c"] integral_sum cdef np.ndarray[float, ndim=2, mode="c"] integral_sqr - if num_type == 1: + + if method == 'norm-coeff': integral_sum = integral.integral_image(image) integral_sqr = integral.integral_image(image**2) # use inversed area for accuracy cdef float inv_area = 1.0 / (template.shape[0] * template.shape[1]) - # calculate template norm according to the following: - # variance ** 2 = 1/K Sigma[(x_k - mean) ** 2] - # = 1/K Sigma[x_k ** 2] - mean ** 2 cdef float template_norm cdef float template_mean = np.mean(template) - if num_type == 0: - template_norm = sqrt((np.std(template) ** 2 + - template_mean ** 2)) / sqrt(inv_area) + 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) @@ -89,18 +94,17 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image, for i in range(result.shape[0] - 1): for j in range(result.shape[1] - 1): num = result[i, j] + i_end = i + template.shape[0] + j_end = j + template.shape[1] + window_mean2 = 0 - if num_type == 1: - t = sum_integral(integral_sum, i, j, - i + template.shape[0], - j + template.shape[1]) + if method == 'norm-coeff': + t = sum_integral(integral_sum, i, j, i_end, j_end) window_mean2 = t * t * inv_area num -= t*template_mean - # calculate squared template window sum in the image - window_sum2 = sum_integral(integral_sqr, i, j, - i + template.shape[0], - j + template.shape[1]) + window_sum2 = sum_integral(integral_sqr, i, j, i_end, j_end) + normed = sqrt(window_sum2 - window_mean2) * template_norm # enforce some limits if fabs(num) < normed: @@ -114,9 +118,7 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image, num = 0 result[i, j] = num # zero boundaries - for i in range(result.shape[0]): - result[i, -1] = 0 - for j in range(result.shape[1]): - result[-1, j] = 0 + result[:, -1] = 0 + result[-1, :] = 0 return result diff --git a/skimage/feature/template.py b/skimage/feature/template.py index b7ff70e1..40306fab 100644 --- a/skimage/feature/template.py +++ b/skimage/feature/template.py @@ -6,9 +6,12 @@ import _template from skimage.util.dtype import _convert -def match_template(image, template, method="norm-coeff"): +def match_template(image, template, method='norm-coeff'): """Finds a template in 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. + Parameters ---------- image : array_like, dtype=float @@ -20,29 +23,26 @@ def match_template(image, template, method="norm-coeff"): T represents the template, I the image and R the result. The summation is done over X = 0..w-1 and Y = 0..h-1 of the template. 'norm-coeff': - R(x, y) = Sigma(X,Y)[T(X, Y).I(x + X, y + Y)] / N - N = sqrt(Sigma(X,Y)[T(X, Y)**2].Sigma(X,Y)[I(x + X, y + Y)**2]) + R(x, y) = Sum(X,Y)[T(X, Y) * I(x + X, y + Y)] / N + N = sqrt(Sum(X,Y)[T(X, Y)**2] * Sum(X,Y)[I(x + X, y + Y)**2]) 'norm-corr': - R(x,y) = Sigma(X,y)[T'(X, Y).I'(x + X, y + Y)] / N - N = sqrt(Sigma(X,y)[T'(X, Y)**2].Sigma(X,Y)[I'(x + X, y + Y)**2]) + R(x,y) = Sum(X,y)[T'(X, Y) * I'(x + X, y + Y)] / N + N = sqrt(Sum(X,y)[T'(X, Y)**2] * Sum(X,Y)[I'(x + X, y + Y)**2]) where: - T'(x, y) = T(X, Y) - 1/(w.h).Sigma(X',Y')[T(X', Y')] + T'(x, y) = T(X, Y) - 1/(w * h) * Sum(X',Y')[T(X', Y')] I'(x + X, y + Y) = I(x + X, y + Y) - - 1/(w.h).Sigma(X',Y')[I(x + X', y + Y')] + - 1/(w * h) * Sum(X',Y')[I(x + X', y + Y')] Returns ------- output : ndarray, dtype=float Correlation results between 0.0 and 1.0, maximum indicating the most probable match. + """ - if method == "norm-corr": - method_num = 0 - elif method == "norm-coeff": - method_num = 1 - else: + 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) - return _template.match_template(image, template, method_num) + return _template.match_template(image, template, method)