From 783c04baec7a8bdcfc6229059ee4877d9bee8a30 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20Sch=C3=B6nberger?= Date: Sun, 8 Dec 2013 18:46:38 +0100 Subject: [PATCH] Fix bug in array padding, template matching --- bento.info | 3 - skimage/_shared/transform.pyx | 1 + skimage/feature/_template.pyx | 97 ----------------- skimage/feature/setup.py | 3 - skimage/feature/template.py | 140 +++++++++++++++++-------- skimage/feature/tests/test_template.py | 3 +- skimage/util/__init__.py | 9 +- skimage/util/arraypad.py | 6 +- 8 files changed, 108 insertions(+), 154 deletions(-) delete mode 100644 skimage/feature/_template.pyx diff --git a/bento.info b/bento.info index 30367057..ff9122fa 100644 --- a/bento.info +++ b/bento.info @@ -39,9 +39,6 @@ Library: Extension: skimage.morphology._pnpoly Sources: skimage/morphology/_pnpoly.pyx - Extension: skimage.feature._template - Sources: - skimage/feature/_template.pyx Extension: skimage.io._plugins._colormixer Sources: skimage/io/_plugins/_colormixer.pyx diff --git a/skimage/_shared/transform.pyx b/skimage/_shared/transform.pyx index 9bdc6824..d77ec583 100644 --- a/skimage/_shared/transform.pyx +++ b/skimage/_shared/transform.pyx @@ -41,4 +41,5 @@ cdef float integrate(float[:, ::1] sat, Py_ssize_t r0, Py_ssize_t c0, if (c0 - 1 >= 0): S -= sat[r1, c0 - 1] + return S diff --git a/skimage/feature/_template.pyx b/skimage/feature/_template.pyx deleted file mode 100644 index 855ece23..00000000 --- a/skimage/feature/_template.pyx +++ /dev/null @@ -1,97 +0,0 @@ -#cython: cdivision=True -#cython: boundscheck=False -#cython: nonecheck=False -#cython: wraparound=False - -""" -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 numpy as np -from scipy.signal import fftconvolve - -cimport numpy as cnp -from libc.math cimport sqrt, fabs -from skimage._shared.transform cimport integrate - - -from skimage.transform import integral - - -def match_template(cnp.ndarray[float, ndim=2, mode="c"] image, - cnp.ndarray[float, ndim=2, mode="c"] template): - - cdef float[:, ::1] corr - cdef float[:, ::1] image_sat - cdef float[:, ::1] image_sqr_sat - cdef float template_mean = np.mean(template) - cdef float template_ssd - cdef float inv_area - cdef Py_ssize_t r, c, r_end, c_end - cdef Py_ssize_t template_rows = template.shape[0] - cdef Py_ssize_t template_cols = template.shape[1] - cdef float den, window_sqr_sum, window_mean_sqr, window_sum - - image_sat = integral.integral_image(image) - image_sqr_sat = 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``. - corr = np.ascontiguousarray(fftconvolve(image, - template[::-1, ::-1], - mode="valid"), - dtype=np.float32) - - - # move window through convolution results, normalizing in the process - for r in range(corr.shape[0]): - for c in range(corr.shape[1]): - # subtract 1 because `i_end` and `c_end` are used for indexing into - # summed-area table, instead of slicing windows of the image. - r_end = r + template_rows - 1 - c_end = c + template_cols - 1 - - window_sum = integrate(image_sat, r, c, r_end, c_end) - window_mean_sqr = window_sum * window_sum * inv_area - window_sqr_sum = integrate(image_sqr_sat, r, c, r_end, c_end) - if window_sqr_sum <= window_mean_sqr: - corr[r, c] = 0 - continue - - den = sqrt((window_sqr_sum - window_mean_sqr) * template_ssd) - corr[r, c] /= den - - return np.asarray(corr) diff --git a/skimage/feature/setup.py b/skimage/feature/setup.py index 7df64c32..915bc351 100644 --- a/skimage/feature/setup.py +++ b/skimage/feature/setup.py @@ -16,7 +16,6 @@ def configuration(parent_package='', top_path=None): cython(['censure_cy.pyx'], working_path=base_path) cython(['_brief_cy.pyx'], working_path=base_path) cython(['_texture.pyx'], working_path=base_path) - cython(['_template.pyx'], working_path=base_path) config.add_extension('corner_cy', sources=['corner_cy.c'], include_dirs=[get_numpy_include_dirs()]) @@ -26,8 +25,6 @@ def configuration(parent_package='', top_path=None): include_dirs=[get_numpy_include_dirs()]) config.add_extension('_texture', sources=['_texture.c'], include_dirs=[get_numpy_include_dirs(), '../_shared']) - config.add_extension('_template', sources=['_template.c'], - include_dirs=[get_numpy_include_dirs(), '../_shared']) return config diff --git a/skimage/feature/template.py b/skimage/feature/template.py index e8f17512..d90fe0c5 100644 --- a/skimage/feature/template.py +++ b/skimage/feature/template.py @@ -1,19 +1,33 @@ -"""template.py - Template matching -""" import numpy as np -from . import _template +from scipy.signal import fftconvolve + +from skimage.util import pad -def match_template(image, template, pad_input=False): +def _window_sum(image, window_shape): + + window_sum = np.cumsum(image, axis=0) + window_sum = (window_sum[window_shape[0]:-1] + - window_sum[:-window_shape[0]-1]) + + window_sum = np.cumsum(window_sum, axis=1) + window_sum = (window_sum[:, window_shape[1]:-1] + - window_sum[:, :-window_shape[1]-1]) + + return window_sum + + +def match_template(image, template, pad_input=False, mode='constant', + constant_values=0): """Match a template to a 2-D image using normalized correlation. The output is an array with values between -1.0 and 1.0, which correspond - to the probability that the template is found at that position. + to the correlation coefficient that the template is found at the position. Parameters ---------- image : array_like - 2-D Image to process. + 2-D Image to process. template : array_like Template to locate. pad_input : bool @@ -22,6 +36,10 @@ def match_template(image, template, pad_input=False): Otherwise, the output is an array with shape `(M - m + 1, N - n + 1)` for an `(M, N)` image and an `(m, n)` template, and matches correspond to origin (top-left corner) of the template. + mode : see `numpy.pad`, optional + Padding mode. + constant_values : see `numpy.pad`, optional + Constant values used in conjunction with ``mode='constant'``. Returns ------- @@ -30,52 +48,92 @@ def match_template(image, template, pad_input=False): `(m, n)` template, the `output` is `(M - m + 1, N - n + 1)` when `pad_input = False` and `(M, N)` when `pad_input = True`. + References + ---------- + .. [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. + Examples -------- >>> template = np.zeros((3, 3)) >>> template[1, 1] = 1 - >>> print(template) - [[ 0. 0. 0.] - [ 0. 1. 0.] - [ 0. 0. 0.]] + >>> template + array([[ 0. 0. 0.] + [ 0. 1. 0.] + [ 0. 0. 0.]]) >>> image = np.zeros((6, 6)) >>> image[1, 1] = 1 >>> image[4, 4] = -1 - >>> print(image) - [[ 0. 0. 0. 0. 0. 0.] - [ 0. 1. 0. 0. 0. 0.] - [ 0. 0. 0. 0. 0. 0.] - [ 0. 0. 0. 0. 0. 0.] - [ 0. 0. 0. 0. -1. 0.] - [ 0. 0. 0. 0. 0. 0.]] + >>> image + array([[ 0. 0. 0. 0. 0. 0.] + [ 0. 1. 0. 0. 0. 0.] + [ 0. 0. 0. 0. 0. 0.] + [ 0. 0. 0. 0. 0. 0.] + [ 0. 0. 0. 0. -1. 0.] + [ 0. 0. 0. 0. 0. 0.]]) >>> result = match_template(image, template) - >>> print(np.round(result, 3)) - [[ 1. -0.125 0. 0. ] - [-0.125 -0.125 0. 0. ] - [ 0. 0. 0.125 0.125] - [ 0. 0. 0.125 -1. ]] + >>> np.round(result, 3) + array([[ 1. -0.125 0. 0. ] + [-0.125 -0.125 0. 0. ] + [ 0. 0. 0.125 0.125] + [ 0. 0. 0.125 -1. ]]) >>> result = match_template(image, template, pad_input=True) - >>> print(np.round(result, 3)) - [[-0.125 -0.125 -0.125 0. 0. 0. ] - [-0.125 1. -0.125 0. 0. 0. ] - [-0.125 -0.125 -0.125 0. 0. 0. ] - [ 0. 0. 0. 0.125 0.125 0.125] - [ 0. 0. 0. 0.125 -1. 0.125] - [ 0. 0. 0. 0.125 0.125 0.125]] + >>> np.round(result, 3) + array([[-0.125 -0.125 -0.125 0. 0. 0. ] + [-0.125 1. -0.125 0. 0. 0. ] + [-0.125 -0.125 -0.125 0. 0. 0. ] + [ 0. 0. 0. 0.125 0.125 0.125] + [ 0. 0. 0. 0.125 -1. 0.125] + [ 0. 0. 0. 0.125 0.125 0.125]]) """ + if np.any(np.less(image.shape, template.shape)): raise ValueError("Image must be larger than template.") - image = np.ascontiguousarray(image, dtype=np.float32) - template = np.ascontiguousarray(template, dtype=np.float32) + + orig_shape = image.shape + + image = np.array(image, dtype=np.float32, copy=False) + + if mode == 'constant': + image = pad(image, pad_width=template.shape, mode=mode, + constant_values=constant_values) + else: + image = pad(image, pad_width=template.shape, mode=mode) + + image_window_sum = _window_sum(image, template.shape) + image_window_sum2 = _window_sum(image**2, template.shape) + + template_area = np.prod(template.shape) + template_ssd = np.sum((template - template.mean())**2) + + xcorr = fftconvolve(image, template[::-1, ::-1], mode="valid")[1:-1, 1:-1] + nom = xcorr - image_window_sum * (template.sum() / template_area) + + denom = image_window_sum2 - image_window_sum**2 / template_area + denom *= template_ssd + np.maximum(denom, 0, out=denom) # sqrt of negative number not allowed + np.sqrt(denom, out=denom) + + response = np.zeros_like(xcorr, dtype=np.float32) + + # avoid zero-division + mask = denom > np.finfo(np.float32).eps + + response[mask] = nom[mask] / denom[mask] if pad_input: - pad_size = tuple(np.array(image.shape) + np.array(template.shape) - 1) - pad_image = np.mean(image) * np.ones(pad_size, dtype=np.float32) - h, w = image.shape - i0, j0 = template.shape - i0 /= 2 - j0 /= 2 - pad_image[i0:i0 + h, j0:j0 + w] = image - image = pad_image - result = _template.match_template(image, template) - return result + r0 = (template.shape[0] - 1) // 2 + r1 = r0 + orig_shape[0] + c0 = (template.shape[1] - 1) // 2 + c1 = c0 + orig_shape[1] + else: + r0 = template.shape[0] - 1 + r1 = r0 + orig_shape[0] - template.shape[0] + 1 + c0 = template.shape[1] - 1 + c1 = c0 + orig_shape[1] - template.shape[1] + 1 + + response = response[r0:r1, c0:c1] + + return response diff --git a/skimage/feature/tests/test_template.py b/skimage/feature/tests/test_template.py index 1b9ff213..e5ce7244 100644 --- a/skimage/feature/tests/test_template.py +++ b/skimage/feature/tests/test_template.py @@ -108,7 +108,8 @@ def test_pad_input(): image[mid, -9:-4] -= template # full min template centered at 12 image[mid, -3:] += template[:, :3] # half max template centered at 18 - result = match_template(image, template, pad_input=True) + result = match_template(image, template, pad_input=True, + constant_values=image.mean()) # get the max and min results. sorted_result = np.argsort(result.flat) diff --git a/skimage/util/__init__.py b/skimage/util/__init__.py index 9cd2bc50..5577e46b 100644 --- a/skimage/util/__init__.py +++ b/skimage/util/__init__.py @@ -3,14 +3,7 @@ from .dtype import (img_as_float, img_as_int, img_as_uint, img_as_ubyte, from .shape import view_as_blocks, view_as_windows from .noise import random_noise -import numpy -ver = numpy.__version__.split('.') -chk = int(ver[0] + ver[1]) -if chk < 18: # Use internal version for numpy versions < 1.8.x - from .arraypad import pad -else: - from numpy import pad -del numpy, ver, chk +from .arraypad import pad from ._regular_grid import regular_grid from .unique import unique_rows diff --git a/skimage/util/arraypad.py b/skimage/util/arraypad.py index 66ad55ea..0bc92d7b 100644 --- a/skimage/util/arraypad.py +++ b/skimage/util/arraypad.py @@ -1027,7 +1027,11 @@ def _normalize_shape(narray, shape): """ normshp = None shapelen = len(np.shape(narray)) - if (isinstance(shape, int)) or shape is None: + + if isinstance(shape, np.ndarray): + shape = shape.tolist() + + if isinstance(shape, (int, float)) or shape is None: normshp = ((shape, shape), ) * shapelen elif (isinstance(shape, (tuple, list)) and isinstance(shape[0], (tuple, list))