From 23311a0993c3aac20449feab7f15e9a35946e861 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sun, 11 Dec 2011 19:16:36 -0500 Subject: [PATCH 01/48] Convert template matching Pull Request to skimage This commit takes the template matching implementation from holtzhau/template (Pull Request #13) and converts the code to use the new package name (scikits.image --> skimage). --- doc/examples/plot_template.py | 70 ++++++++ skimage/detection/__init__.py | 1 + skimage/detection/_template.pyx | 219 +++++++++++++++++++++++ skimage/detection/setup.py | 33 ++++ skimage/detection/template.py | 75 ++++++++ skimage/detection/tests/test_template.py | 46 +++++ skimage/setup.py | 1 + 7 files changed, 445 insertions(+) create mode 100644 doc/examples/plot_template.py create mode 100644 skimage/detection/__init__.py create mode 100644 skimage/detection/_template.pyx create mode 100644 skimage/detection/setup.py create mode 100644 skimage/detection/template.py create mode 100644 skimage/detection/tests/test_template.py diff --git a/doc/examples/plot_template.py b/doc/examples/plot_template.py new file mode 100644 index 00000000..6a98ae0e --- /dev/null +++ b/doc/examples/plot_template.py @@ -0,0 +1,70 @@ +""" +================= +Template Matching +================= + +In this example, we use template matching to identify the occurrence of an +object in an image. The ``match_template`` function uses normalised correlation +techniques to find instances of the "target image" in the "test image". + +The output of ``match_template`` is an image where we can easily identify peaks +by eye. Nevertheless, this example concludes with a simple peak extraction +algorithm to quantify the locations of matches. +""" + +import numpy as np +from skimage.detection import match_template +from numpy.random import randn +import matplotlib.pyplot as plt +import math + +# We first construct a simple image target: +size = 100 +target = np.tri(size) + np.tri(size)[::-1] +target = target.astype(np.float32) + +plt.gray() +plt.imshow(target) +plt.title("Target image") +plt.axis('off') + +# place target in an image at two positions, and add noise. +image = np.zeros((400, 400), dtype=np.float32) +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 + +plt.figure() +plt.imshow(image) +plt.title("Test image") +plt.axis('off') + +# Match the template. +result = match_template(image, target, method='norm-corr') + +plt.figure() +plt.imshow(result) +plt.title("Result from ``match_template``") +plt.axis('off') + +plt.show() + +# peak extraction algorithm. +delta = 5 +found_positions = [] +for i in range(50): + index = np.argmax(result) + y, x = np.unravel_index(index, result.shape) + if not found_positions: + found_positions.append((x, y)) + for position in found_positions: + distance = math.sqrt((x - position[0]) ** 2 + (y - position[1]) ** 2) + if distance > delta: + found_positions.append((x, y)) + result[y, x] = 0 + if len(found_positions) == len(target_positions): + break + +assert np.all(found_positions == target_positions) + diff --git a/skimage/detection/__init__.py b/skimage/detection/__init__.py new file mode 100644 index 00000000..3fdc2389 --- /dev/null +++ b/skimage/detection/__init__.py @@ -0,0 +1 @@ +from template import match_template diff --git a/skimage/detection/_template.pyx b/skimage/detection/_template.pyx new file mode 100644 index 00000000..9c9b9f6f --- /dev/null +++ b/skimage/detection/_template.pyx @@ -0,0 +1,219 @@ +"""template.py - Template matching +""" +import cython +cimport numpy as np +import numpy as np +import cv +from scipy.signal import fftconvolve + +cdef extern from "math.h": + double sqrt(double x) + double fabs(double x) + + +@cython.boundscheck(False) +cdef integral_image(np.ndarray[float, ndim=2, mode="c"] image): + """ + Calculate the summed integral image. + + Parameters + ---------- + image : array_like, dtype=float + Source image. + + Returns + ------- + output : ndarray, dtype=np.double_t + Summed integral image. + """ + cdef np.ndarray[np.double_t, ndim=2, mode="c"] ii = np.zeros((image.shape[0], image.shape[1])) + cdef double s + cdef int x, y + cdef int width, height + height = image.shape[0] + width = image.shape[1] + ii[0, 0] = image[0, 0] + + for y in range(1, height): + ii[y, 0] = image[y, 0] + ii[y - 1, 0] + + for x in range(1, width): + s = 0 + for y in range(0, height): + s += image[y, x] + ii[y, x] = s + ii[y, x - 1] + + return ii + + +@cython.boundscheck(False) +cdef integral_image_sqr(np.ndarray[float, ndim=2, mode="c"] image): + """ + Calculate the squared integral image. + + Parameters + ---------- + image : array_like, dtype=float + Source image. + + Returns + ------- + output : ndarray, dtype=np.double_t + Squared integral image. + """ + cdef np.ndarray[np.double_t, ndim=2, mode="c"] ii2 = np.zeros((image.shape[0], image.shape[1])) + cdef double s + cdef int x, y + cdef int width, height + height = image.shape[0] + width = image.shape[1] + ii2[0, 0] = image[0, 0] * image[0, 0] + + for y in range(1, height): + ii2[y, 0] = image[y, 0] * image[y, 0] + ii2[y - 1, 0] + + for x in range(1, width): + s = 0 + for y in range(0, height): + s += image[y, x] * image[y, x] + ii2[y, x] = s + ii2[y, x - 1] + + return ii2 + + +@cython.boundscheck(False) +cdef integral_images(np.ndarray[float, ndim=2, mode="c"] image): + """ + Calculate the summed and sqared integral image. + + Parameters + ---------- + image : array_like, dtype=float + Source image. + + Returns + ------- + output : tuple (ndarray, ndarray) of type np.double_t + Summed and squared integral image. + """ + cdef np.ndarray[np.double_t, ndim=2, mode="c"] ii = np.zeros((image.shape[0], image.shape[1])) + cdef np.ndarray[np.double_t, ndim=2, mode="c"] ii2 = np.zeros((image.shape[0], image.shape[1])) + cdef double s, s2 + cdef int x, y + cdef int width, height + height = image.shape[0] + width = image.shape[1] + ii[0, 0] = image[0, 0] + ii2[0, 0] = image[0, 0] * image[0, 0] + + for y in range(1, height): + ii[y, 0] = image[y, 0] + ii[y - 1, 0] + ii2[y, 0] = image[y, 0] * image[y, 0] + ii2[y - 1, 0] + + for x in range(1, width): + s = 0 + s2 = 0 + for y in range(0, height): + s += image[y, x] + s2 += image[y, x] * image[y, x] + ii[y, x] = s + ii[y, x - 1] + ii2[y, x] = s2 + ii2[y, x - 1] + + return ii, ii2 + + +@cython.boundscheck(False) +cdef double sum_integral(np.ndarray[np.double_t, ndim=2, mode="c"] sat, + int r0, int c0, int r1, int c1): + """ + Using a summed area table / integral image, calculate the sum + over a given window. + + Parameters + ---------- + sat : ndarray of double_t + Summed area table / integral image. + r0, c0 : int + Top-left corner of block to be summed. + r1, c1 : int + Bottom-right corner of block to be summed. + + Returns + ------- + S : int + Sum over the given window. + """ + cdef double S = 0 + + S += sat[r1, c1] + + if (r0 - 1 >= 0) and (c0 - 1 >= 0): + S += sat[r0 - 1, c0 - 1] + + if (r0 - 1 >= 0): + S -= sat[r0 - 1, c1] + + if (c0 - 1 >= 0): + S -= sat[r1, c0 - 1] + return S + + +@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): + # convolve the image with template by frequency domain multiplication + cdef np.ndarray[np.double_t, ndim=2] result + result = np.ascontiguousarray(fftconvolve(image, np.fliplr(template), mode="valid"), dtype=np.double) + # calculate squared integral images used for normalization + cdef np.ndarray[np.double_t, ndim=2, mode="c"] integral_sum + cdef np.ndarray[np.double_t, ndim=2, mode="c"] integral_sqr + if num_type == 1: + integral_sum, integral_sqr = integral_images(image) + else: + integral_sqr = integral_image_sqr(image) + + # use inversed area for accuracy + cdef double 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 double template_norm + cdef double template_mean = np.mean(template) + + if num_type == 0: + template_norm = sqrt((np.std(template) ** 2 + template_mean ** 2)) / sqrt(inv_area) + else: + template_norm = sqrt((template_mean ** 2)) / sqrt(inv_area) + + # define window of template size in squared integral image + cdef int i, j + cdef double num, window_sum2, window_mean2, normed, t, + # move window through convolution results, normalizing in the process + for i in range(result.shape[0] - 1): + for j in range(result.shape[1] - 1): + num = result[i, j] + window_mean2 = 0 + if num_type == 1: + t = sum_integral(integral_sum, i, j, i + template.shape[0], j + template.shape[1]) + 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]) + normed = sqrt(window_sum2 - window_mean2) * template_norm + # enforce some limits + if fabs(num) < normed: + num /= normed + elif fabs(num) < normed*1.125: + if num > 0: + num = 1 + else: + num = -1 + else: + 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 + return result diff --git a/skimage/detection/setup.py b/skimage/detection/setup.py new file mode 100644 index 00000000..fd297bd4 --- /dev/null +++ b/skimage/detection/setup.py @@ -0,0 +1,33 @@ +#!/usr/bin/env python + +import os +import shutil +import hashlib + +from skimage._build import cython + +base_path = os.path.abspath(os.path.dirname(__file__)) + +def configuration(parent_package='', top_path=None): + from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs + + config = Configuration('transform', parent_package, top_path) + config.add_data_dir('tests') + + cython(['_template.pyx'], working_path=base_path) + + config.add_extension('_template', sources=['_template.c'], + include_dirs=[get_numpy_include_dirs()]) + + return config + +if __name__ == '__main__': + from numpy.distutils.core import setup + setup(maintainer = 'Scikits.Image Developers', + author = 'Scikits.Image Developers', + maintainer_email = 'scikits-image@googlegroups.com', + description = 'Transforms', + url = 'http://stefanv.github.com/scikits.image/', + license = 'SciPy License (BSD Style)', + **(configuration(top_path='').todict()) + ) diff --git a/skimage/detection/template.py b/skimage/detection/template.py new file mode 100644 index 00000000..9ea8dc56 --- /dev/null +++ b/skimage/detection/template.py @@ -0,0 +1,75 @@ +"""template.py - Template matching +""" +import numpy as np +import cv +import _template + +#XXX add to opencv backend once backend system in place +def match_template_cv(image, template, out=None, method="norm-coeff"): + """Finds a template in an image using normalized correlation. + + Parameters + ---------- + image : array_like, dtype=float + Image to process. + template : array_like, dtype=float + Template to locate. + out: array_like, dtype=float, optional + Optional destination. + Returns + ------- + output : ndarray, dtype=float + Correlation results between 0.0 and 1.0, maximum indicating the most probable match. + """ + if out == None: + out = np.empty((image.shape[0] - template.shape[0] + 1,image.shape[1] - template.shape[1] + 1), dtype=image.dtype) + if method == "norm-corr": + cv.MatchTemplate(image, template, out, cv.CV_TM_CCORR_NORMED) + elif method == "norm-corr": + cv.MatchTemplate(image, template, out, cv.CV_TM_CCOEFF_NORMED) + else: + raise ValueError("Unknown template method: %s" % method) + return out + + +def match_template(image, template, method="norm-coeff"): + """Finds a template in an image using normalized correlation. + + Parameters + ---------- + image : array_like, dtype=float + Image to process. + template : array_like, dtype=float + Template to locate. + method: str (default 'norm-coeff') + The correlation method used in scanning. + 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]) + '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]) + where: + T'(x, y) = T(x', y') - 1/(w.h).Sigma(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'')] + + 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: + raise ValueError("Unknown template method: %s" % method) + return _template.match_template(image, template, method_num) + + + + diff --git a/skimage/detection/tests/test_template.py b/skimage/detection/tests/test_template.py new file mode 100644 index 00000000..dfae37cb --- /dev/null +++ b/skimage/detection/tests/test_template.py @@ -0,0 +1,46 @@ +import os.path +import numpy as np +from numpy.testing import * +from skimage import data_dir +from skimage.detection import * +from numpy.random import randn + +def test_template(): + size = 100 + image = np.zeros((400, 400), dtype=np.float32) + target = np.tri(size) + np.tri(size)[::-1] + target = target.astype(np.float32) + 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 + + for method in ["norm-corr", "norm-coeff"]: + result = match_template(image, target, method=method) + delta = 5 + found_positions = [] + # find the targets + for i in range(50): + index = np.argmax(result) + y, x = np.unravel_index(index, result.shape) + if not found_positions: + found_positions.append((x, y)) + for position in found_positions: + distance = np.sqrt((x - position[0]) ** 2 + (y - position[1]) ** 2) + if distance > delta: + found_positions.append((x, y)) + result[y, x] = 0 + if len(found_positions) == len(target_positions): + break + + for x, y in target_positions: + print x, y + found = False + for position in found_positions: + distance = np.sqrt((x - position[0]) ** 2 + (y - position[1]) ** 2) + if distance < delta: + found = True + assert found + +if __name__ == "__main__": + run_module_suite() diff --git a/skimage/setup.py b/skimage/setup.py index c26014f8..752cdcc7 100644 --- a/skimage/setup.py +++ b/skimage/setup.py @@ -16,6 +16,7 @@ def configuration(parent_package='', top_path=None): config.add_subpackage('draw') config.add_subpackage('feature') config.add_subpackage('measure') + config.add_subpackage('detection') def add_test_directories(arg, dirname, fnames): if dirname.split(os.path.sep)[-1] == 'tests': From 424a2b8e5273799318be666e66e87484ba2c5bd9 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sun, 11 Dec 2011 19:54:58 -0500 Subject: [PATCH 02/48] Remove unused imports --- skimage/detection/setup.py | 2 -- skimage/detection/tests/test_template.py | 9 ++++----- 2 files changed, 4 insertions(+), 7 deletions(-) diff --git a/skimage/detection/setup.py b/skimage/detection/setup.py index fd297bd4..52445d0c 100644 --- a/skimage/detection/setup.py +++ b/skimage/detection/setup.py @@ -1,8 +1,6 @@ #!/usr/bin/env python import os -import shutil -import hashlib from skimage._build import cython diff --git a/skimage/detection/tests/test_template.py b/skimage/detection/tests/test_template.py index dfae37cb..c781279b 100644 --- a/skimage/detection/tests/test_template.py +++ b/skimage/detection/tests/test_template.py @@ -1,8 +1,5 @@ -import os.path import numpy as np -from numpy.testing import * -from skimage import data_dir -from skimage.detection import * +from skimage.detection import match_template from numpy.random import randn def test_template(): @@ -43,4 +40,6 @@ def test_template(): assert found if __name__ == "__main__": - run_module_suite() + from numpy import testing + testing.run_module_suite() + From f6b279bff7d867b38468a259687a37c7356d1fb2 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sun, 11 Dec 2011 19:57:40 -0500 Subject: [PATCH 03/48] Fix whitespace --- skimage/detection/_template.pyx | 33 ++++++++++++------------ skimage/detection/setup.py | 1 + skimage/detection/template.py | 14 +++++----- skimage/detection/tests/test_template.py | 6 ++--- 4 files changed, 27 insertions(+), 27 deletions(-) diff --git a/skimage/detection/_template.pyx b/skimage/detection/_template.pyx index 9c9b9f6f..67a462d7 100644 --- a/skimage/detection/_template.pyx +++ b/skimage/detection/_template.pyx @@ -15,12 +15,12 @@ cdef extern from "math.h": cdef integral_image(np.ndarray[float, ndim=2, mode="c"] image): """ Calculate the summed integral image. - + Parameters ---------- image : array_like, dtype=float Source image. - + Returns ------- output : ndarray, dtype=np.double_t @@ -42,7 +42,7 @@ cdef integral_image(np.ndarray[float, ndim=2, mode="c"] image): for y in range(0, height): s += image[y, x] ii[y, x] = s + ii[y, x - 1] - + return ii @@ -50,12 +50,12 @@ cdef integral_image(np.ndarray[float, ndim=2, mode="c"] image): cdef integral_image_sqr(np.ndarray[float, ndim=2, mode="c"] image): """ Calculate the squared integral image. - + Parameters ---------- image : array_like, dtype=float Source image. - + Returns ------- output : ndarray, dtype=np.double_t @@ -77,7 +77,7 @@ cdef integral_image_sqr(np.ndarray[float, ndim=2, mode="c"] image): for y in range(0, height): s += image[y, x] * image[y, x] ii2[y, x] = s + ii2[y, x - 1] - + return ii2 @@ -85,12 +85,12 @@ cdef integral_image_sqr(np.ndarray[float, ndim=2, mode="c"] image): cdef integral_images(np.ndarray[float, ndim=2, mode="c"] image): """ Calculate the summed and sqared integral image. - + Parameters ---------- image : array_like, dtype=float Source image. - + Returns ------- output : tuple (ndarray, ndarray) of type np.double_t @@ -118,12 +118,12 @@ cdef integral_images(np.ndarray[float, ndim=2, mode="c"] image): s2 += image[y, x] * image[y, x] ii[y, x] = s + ii[y, x - 1] ii2[y, x] = s2 + ii2[y, x - 1] - + return ii, ii2 @cython.boundscheck(False) -cdef double sum_integral(np.ndarray[np.double_t, ndim=2, mode="c"] sat, +cdef double sum_integral(np.ndarray[np.double_t, ndim=2, mode="c"] sat, int r0, int c0, int r1, int c1): """ Using a summed area table / integral image, calculate the sum @@ -178,15 +178,15 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image, # variance ** 2 = 1/K Sigma[(x_k - mean) ** 2] = 1/K Sigma[x_k ** 2] - mean ** 2 cdef double template_norm cdef double template_mean = np.mean(template) - + if num_type == 0: template_norm = sqrt((np.std(template) ** 2 + template_mean ** 2)) / sqrt(inv_area) else: template_norm = sqrt((template_mean ** 2)) / sqrt(inv_area) - + # define window of template size in squared integral image cdef int i, j - cdef double num, window_sum2, window_mean2, normed, t, + cdef double num, window_sum2, window_mean2, normed, t, # move window through convolution results, normalizing in the process for i in range(result.shape[0] - 1): for j in range(result.shape[1] - 1): @@ -196,7 +196,7 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image, t = sum_integral(integral_sum, i, j, i + template.shape[0], j + template.shape[1]) 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]) normed = sqrt(window_sum2 - window_mean2) * template_norm @@ -207,7 +207,7 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image, if num > 0: num = 1 else: - num = -1 + num = -1 else: num = 0 result[i, j] = num @@ -215,5 +215,6 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image, for i in range(result.shape[0]): result[i, -1] = 0 for j in range(result.shape[1]): - result[-1, j] = 0 + result[-1, j] = 0 return result + diff --git a/skimage/detection/setup.py b/skimage/detection/setup.py index 52445d0c..6295cba5 100644 --- a/skimage/detection/setup.py +++ b/skimage/detection/setup.py @@ -29,3 +29,4 @@ if __name__ == '__main__': license = 'SciPy License (BSD Style)', **(configuration(top_path='').todict()) ) + diff --git a/skimage/detection/template.py b/skimage/detection/template.py index 9ea8dc56..eb780799 100644 --- a/skimage/detection/template.py +++ b/skimage/detection/template.py @@ -4,6 +4,7 @@ import numpy as np import cv import _template + #XXX add to opencv backend once backend system in place def match_template_cv(image, template, out=None, method="norm-coeff"): """Finds a template in an image using normalized correlation. @@ -43,17 +44,17 @@ def match_template(image, template, method="norm-coeff"): Template to locate. method: str (default 'norm-coeff') The correlation method used in scanning. - T represents the template, I the image and R the result. + 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': + '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]) - 'norm-corr': - R(x,y) = Sigma(x',y)[T'(x', y').I'(x + x', y + y')] / N + '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]) where: T'(x, y) = T(x', y') - 1/(w.h).Sigma(x'',y'')[T(x'', y'')] - I'(x + x', y + y') = I(x + x', y + y') - + I'(x + x', y + y') = I(x + x', y + y') - 1/(w.h).Sigma(x'',y'')[I(x + x'', y + y'')] Returns @@ -70,6 +71,3 @@ def match_template(image, template, method="norm-coeff"): raise ValueError("Unknown template method: %s" % method) return _template.match_template(image, template, method_num) - - - diff --git a/skimage/detection/tests/test_template.py b/skimage/detection/tests/test_template.py index c781279b..987a3496 100644 --- a/skimage/detection/tests/test_template.py +++ b/skimage/detection/tests/test_template.py @@ -11,14 +11,14 @@ def test_template(): for x, y in target_positions: image[x:x+size, y:y+size] = target image += randn(400, 400)*2 - + for method in ["norm-corr", "norm-coeff"]: result = match_template(image, target, method=method) delta = 5 found_positions = [] # find the targets for i in range(50): - index = np.argmax(result) + index = np.argmax(result) y, x = np.unravel_index(index, result.shape) if not found_positions: found_positions.append((x, y)) @@ -38,7 +38,7 @@ def test_template(): if distance < delta: found = True assert found - + if __name__ == "__main__": from numpy import testing testing.run_module_suite() From b79a8dd437df992aa3275339dc20c2a4fae9303f Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sun, 11 Dec 2011 20:04:08 -0500 Subject: [PATCH 04/48] Update setup file to new repo and fix typos --- skimage/detection/setup.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/skimage/detection/setup.py b/skimage/detection/setup.py index 6295cba5..a5b7a6b4 100644 --- a/skimage/detection/setup.py +++ b/skimage/detection/setup.py @@ -1,7 +1,6 @@ #!/usr/bin/env python import os - from skimage._build import cython base_path = os.path.abspath(os.path.dirname(__file__)) @@ -9,7 +8,7 @@ base_path = os.path.abspath(os.path.dirname(__file__)) def configuration(parent_package='', top_path=None): from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs - config = Configuration('transform', parent_package, top_path) + config = Configuration('detection', parent_package, top_path) config.add_data_dir('tests') cython(['_template.pyx'], working_path=base_path) @@ -21,11 +20,11 @@ def configuration(parent_package='', top_path=None): if __name__ == '__main__': from numpy.distutils.core import setup - setup(maintainer = 'Scikits.Image Developers', - author = 'Scikits.Image Developers', + setup(maintainer = 'scikits-image Developers', + author = 'scikits-image Developers', maintainer_email = 'scikits-image@googlegroups.com', - description = 'Transforms', - url = 'http://stefanv.github.com/scikits.image/', + description = 'detection', + url = 'https://github.com/scikits-image/scikits-image', license = 'SciPy License (BSD Style)', **(configuration(top_path='').todict()) ) From f8e0478542419c9266249367c34aa44d70976666 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sun, 11 Dec 2011 20:05:02 -0500 Subject: [PATCH 05/48] Remove dependency on math module --- doc/examples/plot_template.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/doc/examples/plot_template.py b/doc/examples/plot_template.py index 6a98ae0e..28bfbcfe 100644 --- a/doc/examples/plot_template.py +++ b/doc/examples/plot_template.py @@ -16,7 +16,6 @@ import numpy as np from skimage.detection import match_template from numpy.random import randn import matplotlib.pyplot as plt -import math # We first construct a simple image target: size = 100 @@ -59,7 +58,7 @@ for i in range(50): if not found_positions: found_positions.append((x, y)) for position in found_positions: - distance = math.sqrt((x - position[0]) ** 2 + (y - position[1]) ** 2) + distance = np.sqrt((x - position[0]) ** 2 + (y - position[1]) ** 2) if distance > delta: found_positions.append((x, y)) result[y, x] = 0 From ae85fa8b153181023c9c67b650927cf35600415f Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sun, 11 Dec 2011 20:15:31 -0500 Subject: [PATCH 06/48] Wrap long lines --- skimage/detection/_template.pyx | 11 ++++++++--- skimage/detection/template.py | 7 +++++-- skimage/detection/tests/test_template.py | 6 ++++-- 3 files changed, 17 insertions(+), 7 deletions(-) diff --git a/skimage/detection/_template.pyx b/skimage/detection/_template.pyx index 67a462d7..90bfe666 100644 --- a/skimage/detection/_template.pyx +++ b/skimage/detection/_template.pyx @@ -163,7 +163,8 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image, np.ndarray[float, ndim=2, mode="c"] template, int num_type): # convolve the image with template by frequency domain multiplication cdef np.ndarray[np.double_t, ndim=2] result - result = np.ascontiguousarray(fftconvolve(image, np.fliplr(template), mode="valid"), dtype=np.double) + result = np.ascontiguousarray(fftconvolve(image, np.fliplr(template), + mode="valid"), dtype=np.double) # calculate squared integral images used for normalization cdef np.ndarray[np.double_t, ndim=2, mode="c"] integral_sum cdef np.ndarray[np.double_t, ndim=2, mode="c"] integral_sqr @@ -193,12 +194,16 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image, num = result[i, j] window_mean2 = 0 if num_type == 1: - t = sum_integral(integral_sum, i, j, i + template.shape[0], j + template.shape[1]) + t = sum_integral(integral_sum, i, j, + i + template.shape[0], + j + template.shape[1]) 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 + template.shape[0], + j + template.shape[1]) normed = sqrt(window_sum2 - window_mean2) * template_norm # enforce some limits if fabs(num) < normed: diff --git a/skimage/detection/template.py b/skimage/detection/template.py index eb780799..3d8b9a0c 100644 --- a/skimage/detection/template.py +++ b/skimage/detection/template.py @@ -20,10 +20,13 @@ def match_template_cv(image, template, out=None, method="norm-coeff"): Returns ------- output : ndarray, dtype=float - Correlation results between 0.0 and 1.0, maximum indicating the most probable match. + Correlation results between 0.0 and 1.0, maximum indicating the most + probable match. """ if out == None: - out = np.empty((image.shape[0] - template.shape[0] + 1,image.shape[1] - template.shape[1] + 1), dtype=image.dtype) + out = np.empty((image.shape[0] - template.shape[0] + 1, + image.shape[1] - template.shape[1] + 1), + dtype=image.dtype) if method == "norm-corr": cv.MatchTemplate(image, template, out, cv.CV_TM_CCORR_NORMED) elif method == "norm-corr": diff --git a/skimage/detection/tests/test_template.py b/skimage/detection/tests/test_template.py index 987a3496..b5cc9bcd 100644 --- a/skimage/detection/tests/test_template.py +++ b/skimage/detection/tests/test_template.py @@ -23,7 +23,8 @@ def test_template(): if not found_positions: found_positions.append((x, y)) for position in found_positions: - distance = np.sqrt((x - position[0]) ** 2 + (y - position[1]) ** 2) + distance = np.sqrt((x - position[0]) ** 2 + + (y - position[1]) ** 2) if distance > delta: found_positions.append((x, y)) result[y, x] = 0 @@ -34,7 +35,8 @@ def test_template(): print x, y found = False for position in found_positions: - distance = np.sqrt((x - position[0]) ** 2 + (y - position[1]) ** 2) + distance = np.sqrt((x - position[0]) ** 2 + + (y - position[1]) ** 2) if distance < delta: found = True assert found From 19f1021f60c32ffdb8281768309c7801b604b568 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sun, 11 Dec 2011 22:19:56 -0500 Subject: [PATCH 07/48] Clean up docstrings --- skimage/detection/_template.pyx | 2 +- skimage/detection/template.py | 16 ++++++++-------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/skimage/detection/_template.pyx b/skimage/detection/_template.pyx index 90bfe666..fb9d4028 100644 --- a/skimage/detection/_template.pyx +++ b/skimage/detection/_template.pyx @@ -84,7 +84,7 @@ cdef integral_image_sqr(np.ndarray[float, ndim=2, mode="c"] image): @cython.boundscheck(False) cdef integral_images(np.ndarray[float, ndim=2, mode="c"] image): """ - Calculate the summed and sqared integral image. + Calculate the summed and squared integral image. Parameters ---------- diff --git a/skimage/detection/template.py b/skimage/detection/template.py index 3d8b9a0c..91412e94 100644 --- a/skimage/detection/template.py +++ b/skimage/detection/template.py @@ -48,17 +48,17 @@ def match_template(image, template, method="norm-coeff"): method: str (default 'norm-coeff') The correlation method used in scanning. 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. + 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) = 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]) '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) = 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]) where: - T'(x, y) = T(x', y') - 1/(w.h).Sigma(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'')] + T'(x, y) = T(X, Y) - 1/(w.h).Sigma(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')] Returns ------- From 49349a07c6ca1dd94f157793996597b20704c244 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sun, 11 Dec 2011 22:20:52 -0500 Subject: [PATCH 08/48] Fix copy-paste bug --- skimage/detection/template.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/detection/template.py b/skimage/detection/template.py index 91412e94..e0ad01ef 100644 --- a/skimage/detection/template.py +++ b/skimage/detection/template.py @@ -29,7 +29,7 @@ def match_template_cv(image, template, out=None, method="norm-coeff"): dtype=image.dtype) if method == "norm-corr": cv.MatchTemplate(image, template, out, cv.CV_TM_CCORR_NORMED) - elif method == "norm-corr": + elif method == "norm-coeff": cv.MatchTemplate(image, template, out, cv.CV_TM_CCOEFF_NORMED) else: raise ValueError("Unknown template method: %s" % method) From d43049da713df6bab2e56a755f9f50f764f5109b Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sun, 11 Dec 2011 23:05:57 -0500 Subject: [PATCH 09/48] Wrap long lines --- skimage/detection/_template.pyx | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/skimage/detection/_template.pyx b/skimage/detection/_template.pyx index fb9d4028..3c6a726f 100644 --- a/skimage/detection/_template.pyx +++ b/skimage/detection/_template.pyx @@ -176,12 +176,14 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image, # use inversed area for accuracy cdef double 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 + # variance ** 2 = 1/K Sigma[(x_k - mean) ** 2] + # = 1/K Sigma[x_k ** 2] - mean ** 2 cdef double template_norm cdef double template_mean = np.mean(template) if num_type == 0: - template_norm = sqrt((np.std(template) ** 2 + template_mean ** 2)) / sqrt(inv_area) + template_norm = sqrt((np.std(template) ** 2 + + template_mean ** 2)) / sqrt(inv_area) else: template_norm = sqrt((template_mean ** 2)) / sqrt(inv_area) From 4ab4213749cb90cec574c8e066b921e6e162279f Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Mon, 12 Dec 2011 20:19:15 -0500 Subject: [PATCH 10/48] Change OpenCV import to make it an optional dependency --- skimage/detection/_template.pyx | 2 +- skimage/detection/template.py | 10 +++++++++- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/skimage/detection/_template.pyx b/skimage/detection/_template.pyx index 3c6a726f..eb5a9733 100644 --- a/skimage/detection/_template.pyx +++ b/skimage/detection/_template.pyx @@ -3,9 +3,9 @@ import cython cimport numpy as np import numpy as np -import cv from scipy.signal import fftconvolve + cdef extern from "math.h": double sqrt(double x) double fabs(double x) diff --git a/skimage/detection/template.py b/skimage/detection/template.py index e0ad01ef..6bf02896 100644 --- a/skimage/detection/template.py +++ b/skimage/detection/template.py @@ -1,9 +1,15 @@ """template.py - Template matching """ import numpy as np -import cv import _template +try: + import cv + opencv_available = True +except ImportError: + opencv_available = False + + #XXX add to opencv backend once backend system in place def match_template_cv(image, template, out=None, method="norm-coeff"): @@ -23,6 +29,8 @@ def match_template_cv(image, template, out=None, method="norm-coeff"): Correlation results between 0.0 and 1.0, maximum indicating the most probable match. """ + if not opencv_available: + raise ImportError("Opencv 2.0+ required") if out == None: out = np.empty((image.shape[0] - template.shape[0] + 1, image.shape[1] - template.shape[1] + 1), From 56611198e486bde21e4eff954b8c16f00713d8fc Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Mon, 12 Dec 2011 20:25:55 -0500 Subject: [PATCH 11/48] PEP8: add spacing around operators --- skimage/detection/tests/test_template.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/skimage/detection/tests/test_template.py b/skimage/detection/tests/test_template.py index b5cc9bcd..713b6e96 100644 --- a/skimage/detection/tests/test_template.py +++ b/skimage/detection/tests/test_template.py @@ -9,8 +9,8 @@ def test_template(): target = target.astype(np.float32) 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 + image[x:x + size, y:y + size] = target + image += randn(400, 400) * 2 for method in ["norm-corr", "norm-coeff"]: result = match_template(image, target, method=method) From 2b368aecb105d3421aa24220e84042aa1da8b9d7 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Mon, 12 Dec 2011 20:32:20 -0500 Subject: [PATCH 12/48] Remove unused `integral_image` function in _template.pyx --- skimage/detection/_template.pyx | 35 --------------------------------- 1 file changed, 35 deletions(-) diff --git a/skimage/detection/_template.pyx b/skimage/detection/_template.pyx index eb5a9733..3300bf88 100644 --- a/skimage/detection/_template.pyx +++ b/skimage/detection/_template.pyx @@ -11,41 +11,6 @@ cdef extern from "math.h": double fabs(double x) -@cython.boundscheck(False) -cdef integral_image(np.ndarray[float, ndim=2, mode="c"] image): - """ - Calculate the summed integral image. - - Parameters - ---------- - image : array_like, dtype=float - Source image. - - Returns - ------- - output : ndarray, dtype=np.double_t - Summed integral image. - """ - cdef np.ndarray[np.double_t, ndim=2, mode="c"] ii = np.zeros((image.shape[0], image.shape[1])) - cdef double s - cdef int x, y - cdef int width, height - height = image.shape[0] - width = image.shape[1] - ii[0, 0] = image[0, 0] - - for y in range(1, height): - ii[y, 0] = image[y, 0] + ii[y - 1, 0] - - for x in range(1, width): - s = 0 - for y in range(0, height): - s += image[y, x] - ii[y, x] = s + ii[y, x - 1] - - return ii - - @cython.boundscheck(False) cdef integral_image_sqr(np.ndarray[float, ndim=2, mode="c"] image): """ From 6272c312c416a53cff9caa8fb9c25d3f17f5618e Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Mon, 12 Dec 2011 23:02:44 -0500 Subject: [PATCH 13/48] Simplify _template.pyx using integral_image from transform subpackage. Remove `integral_images` and `integral_image_sqr` from _template.pyx in favor of calls to `skimage.transform.integral_image`. This change required `match_template` arguments ("image" and "template") to be changed from float to double. After this change, the template test runs about 25% slower. --- doc/examples/plot_template.py | 3 +- skimage/detection/_template.pyx | 87 ++---------------------- skimage/detection/tests/test_template.py | 4 +- 3 files changed, 9 insertions(+), 85 deletions(-) diff --git a/doc/examples/plot_template.py b/doc/examples/plot_template.py index 28bfbcfe..7894cdbc 100644 --- a/doc/examples/plot_template.py +++ b/doc/examples/plot_template.py @@ -20,7 +20,6 @@ import matplotlib.pyplot as plt # We first construct a simple image target: size = 100 target = np.tri(size) + np.tri(size)[::-1] -target = target.astype(np.float32) plt.gray() plt.imshow(target) @@ -28,7 +27,7 @@ plt.title("Target image") plt.axis('off') # place target in an image at two positions, and add noise. -image = np.zeros((400, 400), dtype=np.float32) +image = np.zeros((400, 400)) target_positions = [(50, 50), (200, 200)] for x, y in target_positions: image[x:x+size, y:y+size] = target diff --git a/skimage/detection/_template.pyx b/skimage/detection/_template.pyx index 3300bf88..40ed261b 100644 --- a/skimage/detection/_template.pyx +++ b/skimage/detection/_template.pyx @@ -4,6 +4,7 @@ import cython cimport numpy as np import numpy as np from scipy.signal import fftconvolve +from skimage.transform import integral cdef extern from "math.h": @@ -11,82 +12,6 @@ cdef extern from "math.h": double fabs(double x) -@cython.boundscheck(False) -cdef integral_image_sqr(np.ndarray[float, ndim=2, mode="c"] image): - """ - Calculate the squared integral image. - - Parameters - ---------- - image : array_like, dtype=float - Source image. - - Returns - ------- - output : ndarray, dtype=np.double_t - Squared integral image. - """ - cdef np.ndarray[np.double_t, ndim=2, mode="c"] ii2 = np.zeros((image.shape[0], image.shape[1])) - cdef double s - cdef int x, y - cdef int width, height - height = image.shape[0] - width = image.shape[1] - ii2[0, 0] = image[0, 0] * image[0, 0] - - for y in range(1, height): - ii2[y, 0] = image[y, 0] * image[y, 0] + ii2[y - 1, 0] - - for x in range(1, width): - s = 0 - for y in range(0, height): - s += image[y, x] * image[y, x] - ii2[y, x] = s + ii2[y, x - 1] - - return ii2 - - -@cython.boundscheck(False) -cdef integral_images(np.ndarray[float, ndim=2, mode="c"] image): - """ - Calculate the summed and squared integral image. - - Parameters - ---------- - image : array_like, dtype=float - Source image. - - Returns - ------- - output : tuple (ndarray, ndarray) of type np.double_t - Summed and squared integral image. - """ - cdef np.ndarray[np.double_t, ndim=2, mode="c"] ii = np.zeros((image.shape[0], image.shape[1])) - cdef np.ndarray[np.double_t, ndim=2, mode="c"] ii2 = np.zeros((image.shape[0], image.shape[1])) - cdef double s, s2 - cdef int x, y - cdef int width, height - height = image.shape[0] - width = image.shape[1] - ii[0, 0] = image[0, 0] - ii2[0, 0] = image[0, 0] * image[0, 0] - - for y in range(1, height): - ii[y, 0] = image[y, 0] + ii[y - 1, 0] - ii2[y, 0] = image[y, 0] * image[y, 0] + ii2[y - 1, 0] - - for x in range(1, width): - s = 0 - s2 = 0 - for y in range(0, height): - s += image[y, x] - s2 += image[y, x] * image[y, x] - ii[y, x] = s + ii[y, x - 1] - ii2[y, x] = s2 + ii2[y, x - 1] - - return ii, ii2 - - @cython.boundscheck(False) cdef double sum_integral(np.ndarray[np.double_t, ndim=2, mode="c"] sat, int r0, int c0, int r1, int c1): @@ -124,8 +49,9 @@ cdef double sum_integral(np.ndarray[np.double_t, 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): +def match_template(np.ndarray[np.double_t, ndim=2, mode="c"] image, + np.ndarray[np.double_t, ndim=2, mode="c"] template, + int num_type): # convolve the image with template by frequency domain multiplication cdef np.ndarray[np.double_t, ndim=2] result result = np.ascontiguousarray(fftconvolve(image, np.fliplr(template), @@ -134,9 +60,8 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image, cdef np.ndarray[np.double_t, ndim=2, mode="c"] integral_sum cdef np.ndarray[np.double_t, ndim=2, mode="c"] integral_sqr if num_type == 1: - integral_sum, integral_sqr = integral_images(image) - else: - integral_sqr = integral_image_sqr(image) + integral_sum = integral.integral_image(image) + integral_sqr = integral.integral_image(image**2) # use inversed area for accuracy cdef double inv_area = 1.0 / (template.shape[0] * template.shape[1]) diff --git a/skimage/detection/tests/test_template.py b/skimage/detection/tests/test_template.py index 713b6e96..a3dfbcde 100644 --- a/skimage/detection/tests/test_template.py +++ b/skimage/detection/tests/test_template.py @@ -2,11 +2,11 @@ import numpy as np from skimage.detection import match_template from numpy.random import randn + def test_template(): size = 100 - image = np.zeros((400, 400), dtype=np.float32) + image = np.zeros((400, 400)) target = np.tri(size) + np.tri(size)[::-1] - target = target.astype(np.float32) target_positions = [(50, 50), (200, 200)] for x, y in target_positions: image[x:x + size, y:y + size] = target From acc35c843b0afec51cd98ea85ce96836be61a990 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Mon, 12 Dec 2011 23:22:45 -0500 Subject: [PATCH 14/48] Add comment about code duplication. --- skimage/detection/_template.pyx | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/skimage/detection/_template.pyx b/skimage/detection/_template.pyx index 40ed261b..f3e89f8a 100644 --- a/skimage/detection/_template.pyx +++ b/skimage/detection/_template.pyx @@ -19,6 +19,10 @@ cdef double sum_integral(np.ndarray[np.double_t, ndim=2, mode="c"] sat, Using a summed area table / integral image, calculate the sum over a given window. + This function is the same as the `integrate` function in + `skimage.transform.integrate`, but this Cython version significantly + speeds up the code. + Parameters ---------- sat : ndarray of double_t From 12b39dae5c6af9f923ee824fd311f0d5d6849f02 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sun, 18 Dec 2011 12:27:33 -0500 Subject: [PATCH 15/48] Fix assert failure in example --- doc/examples/plot_template.py | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/examples/plot_template.py b/doc/examples/plot_template.py index 7894cdbc..861c6225 100644 --- a/doc/examples/plot_template.py +++ b/doc/examples/plot_template.py @@ -64,5 +64,6 @@ for i in range(50): if len(found_positions) == len(target_positions): break +found_positions = np.sort(found_positions) assert np.all(found_positions == target_positions) From e6098e140b011f26acad1960e2c3c70c338a2a6c Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sun, 18 Dec 2011 12:40:49 -0500 Subject: [PATCH 16/48] Replace assert statement with plot to show matches --- doc/examples/plot_template.py | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/doc/examples/plot_template.py b/doc/examples/plot_template.py index 861c6225..8de50fd6 100644 --- a/doc/examples/plot_template.py +++ b/doc/examples/plot_template.py @@ -9,7 +9,7 @@ techniques to find instances of the "target image" in the "test image". The output of ``match_template`` is an image where we can easily identify peaks by eye. Nevertheless, this example concludes with a simple peak extraction -algorithm to quantify the locations of matches. +algorithm to quantify the locations of matches (marked in red). """ import numpy as np @@ -21,7 +21,10 @@ import matplotlib.pyplot as plt size = 100 target = np.tri(size) + np.tri(size)[::-1] -plt.gray() +#plt.gray() +plt.figure(figsize=(9, 3)) + +plt.subplot(1, 3, 1) plt.imshow(target) plt.title("Target image") plt.axis('off') @@ -33,7 +36,7 @@ for x, y in target_positions: image[x:x+size, y:y+size] = target image += randn(400, 400)*2 -plt.figure() +plt.subplot(1, 3, 2) plt.imshow(image) plt.title("Test image") plt.axis('off') @@ -41,13 +44,11 @@ plt.axis('off') # Match the template. result = match_template(image, target, method='norm-corr') -plt.figure() +plt.subplot(1, 3, 3) plt.imshow(result) -plt.title("Result from ``match_template``") +plt.title("Result from\n``match_template``") plt.axis('off') -plt.show() - # peak extraction algorithm. delta = 5 found_positions = [] @@ -64,6 +65,7 @@ for i in range(50): if len(found_positions) == len(target_positions): break -found_positions = np.sort(found_positions) -assert np.all(found_positions == target_positions) - +x_found, y_found = np.transpose(found_positions) +plt.plot(x_found, y_found, 'ro') +plt.autoscale(tight=True) +plt.show() From 01d66fc5012431b2f40b66763606b95b530b7cc9 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sun, 18 Dec 2011 12:52:40 -0500 Subject: [PATCH 17/48] Reorganize example so that all plotting code is at the end --- doc/examples/plot_template.py | 42 +++++++++++++++++------------------ 1 file changed, 20 insertions(+), 22 deletions(-) diff --git a/doc/examples/plot_template.py b/doc/examples/plot_template.py index 8de50fd6..e4a51787 100644 --- a/doc/examples/plot_template.py +++ b/doc/examples/plot_template.py @@ -8,8 +8,8 @@ object in an image. The ``match_template`` function uses normalised correlation techniques to find instances of the "target image" in the "test image". The output of ``match_template`` is an image where we can easily identify peaks -by eye. Nevertheless, this example concludes with a simple peak extraction -algorithm to quantify the locations of matches (marked in red). +by eye. We mark the locations of matches (red dots), which are detected using +a simple peak extraction algorithm. """ import numpy as np @@ -20,15 +20,6 @@ import matplotlib.pyplot as plt # We first construct a simple image target: size = 100 target = np.tri(size) + np.tri(size)[::-1] - -#plt.gray() -plt.figure(figsize=(9, 3)) - -plt.subplot(1, 3, 1) -plt.imshow(target) -plt.title("Target image") -plt.axis('off') - # place target in an image at two positions, and add noise. image = np.zeros((400, 400)) target_positions = [(50, 50), (200, 200)] @@ -36,19 +27,9 @@ for x, y in target_positions: image[x:x+size, y:y+size] = target image += randn(400, 400)*2 -plt.subplot(1, 3, 2) -plt.imshow(image) -plt.title("Test image") -plt.axis('off') - # Match the template. result = match_template(image, target, method='norm-corr') -plt.subplot(1, 3, 3) -plt.imshow(result) -plt.title("Result from\n``match_template``") -plt.axis('off') - # peak extraction algorithm. delta = 5 found_positions = [] @@ -64,8 +45,25 @@ for i in range(50): result[y, x] = 0 if len(found_positions) == len(target_positions): break - x_found, y_found = np.transpose(found_positions) + +plt.gray() + +plt.subplot(1, 3, 1) +plt.imshow(target) +plt.title("Target image") +plt.axis('off') + +plt.subplot(1, 3, 2) +plt.imshow(image) +plt.title("Test image") +plt.axis('off') + +plt.subplot(1, 3, 3) +plt.imshow(result) plt.plot(x_found, y_found, 'ro') +plt.title("Result from\n``match_template``") plt.autoscale(tight=True) +plt.axis('off') + plt.show() From e8461e22dd080fd48d24ca9c2eb0479a6a0298f6 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sun, 18 Dec 2011 13:37:31 -0500 Subject: [PATCH 18/48] Move template matching to feature subpackage --- doc/examples/plot_template.py | 2 +- skimage/detection/__init__.py | 1 - skimage/detection/setup.py | 31 ------------------- skimage/feature/__init__.py | 1 + skimage/{detection => feature}/_template.pyx | 0 skimage/feature/setup.py | 3 ++ skimage/{detection => feature}/template.py | 0 .../tests/test_template.py | 2 +- skimage/setup.py | 1 - 9 files changed, 6 insertions(+), 35 deletions(-) delete mode 100644 skimage/detection/__init__.py delete mode 100644 skimage/detection/setup.py rename skimage/{detection => feature}/_template.pyx (100%) rename skimage/{detection => feature}/template.py (100%) rename skimage/{detection => feature}/tests/test_template.py (97%) diff --git a/doc/examples/plot_template.py b/doc/examples/plot_template.py index e4a51787..3690b18e 100644 --- a/doc/examples/plot_template.py +++ b/doc/examples/plot_template.py @@ -13,7 +13,7 @@ a simple peak extraction algorithm. """ import numpy as np -from skimage.detection import match_template +from skimage.feature import match_template from numpy.random import randn import matplotlib.pyplot as plt diff --git a/skimage/detection/__init__.py b/skimage/detection/__init__.py deleted file mode 100644 index 3fdc2389..00000000 --- a/skimage/detection/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from template import match_template diff --git a/skimage/detection/setup.py b/skimage/detection/setup.py deleted file mode 100644 index a5b7a6b4..00000000 --- a/skimage/detection/setup.py +++ /dev/null @@ -1,31 +0,0 @@ -#!/usr/bin/env python - -import os -from skimage._build import cython - -base_path = os.path.abspath(os.path.dirname(__file__)) - -def configuration(parent_package='', top_path=None): - from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs - - config = Configuration('detection', parent_package, top_path) - config.add_data_dir('tests') - - cython(['_template.pyx'], working_path=base_path) - - config.add_extension('_template', sources=['_template.c'], - include_dirs=[get_numpy_include_dirs()]) - - return config - -if __name__ == '__main__': - from numpy.distutils.core import setup - setup(maintainer = 'scikits-image Developers', - author = 'scikits-image Developers', - maintainer_email = 'scikits-image@googlegroups.com', - description = 'detection', - url = 'https://github.com/scikits-image/scikits-image', - license = 'SciPy License (BSD Style)', - **(configuration(top_path='').todict()) - ) - diff --git a/skimage/feature/__init__.py b/skimage/feature/__init__.py index 4e5d6324..70c154b8 100644 --- a/skimage/feature/__init__.py +++ b/skimage/feature/__init__.py @@ -2,3 +2,4 @@ from .hog import hog from .greycomatrix import greycomatrix, greycoprops from .peak import peak_local_max from .harris import harris +from .template import match_template diff --git a/skimage/detection/_template.pyx b/skimage/feature/_template.pyx similarity index 100% rename from skimage/detection/_template.pyx rename to skimage/feature/_template.pyx diff --git a/skimage/feature/setup.py b/skimage/feature/setup.py index 626b2f5a..13d4fae5 100644 --- a/skimage/feature/setup.py +++ b/skimage/feature/setup.py @@ -12,9 +12,12 @@ def configuration(parent_package='', top_path=None): config.add_data_dir('tests') cython(['_greycomatrix.pyx'], working_path=base_path) + cython(['_template.pyx'], working_path=base_path) config.add_extension('_greycomatrix', sources=['_greycomatrix.c'], include_dirs=[get_numpy_include_dirs()]) + config.add_extension('_template', sources=['_template.c'], + include_dirs=[get_numpy_include_dirs()]) return config diff --git a/skimage/detection/template.py b/skimage/feature/template.py similarity index 100% rename from skimage/detection/template.py rename to skimage/feature/template.py diff --git a/skimage/detection/tests/test_template.py b/skimage/feature/tests/test_template.py similarity index 97% rename from skimage/detection/tests/test_template.py rename to skimage/feature/tests/test_template.py index a3dfbcde..e92672da 100644 --- a/skimage/detection/tests/test_template.py +++ b/skimage/feature/tests/test_template.py @@ -1,5 +1,5 @@ import numpy as np -from skimage.detection import match_template +from skimage.feature import match_template from numpy.random import randn diff --git a/skimage/setup.py b/skimage/setup.py index 752cdcc7..c26014f8 100644 --- a/skimage/setup.py +++ b/skimage/setup.py @@ -16,7 +16,6 @@ def configuration(parent_package='', top_path=None): config.add_subpackage('draw') config.add_subpackage('feature') config.add_subpackage('measure') - config.add_subpackage('detection') def add_test_directories(arg, dirname, fnames): if dirname.split(os.path.sep)[-1] == 'tests': From 2a00e24b12698b1fb6786cfba40b9bb2f9faebbf Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Mon, 26 Dec 2011 11:18:41 -0800 Subject: [PATCH 19/48] Change _template.match_template variables to float. Template and image in test are converted to float32 before passing to `match_template`. This change is temporary: `match_template` should convert these variables internally. --- skimage/feature/_template.pyx | 31 +++++++++++++------------- skimage/feature/tests/test_template.py | 3 ++- 2 files changed, 18 insertions(+), 16 deletions(-) diff --git a/skimage/feature/_template.pyx b/skimage/feature/_template.pyx index f3e89f8a..d64edf04 100644 --- a/skimage/feature/_template.pyx +++ b/skimage/feature/_template.pyx @@ -8,12 +8,12 @@ from skimage.transform import integral cdef extern from "math.h": - double sqrt(double x) - double fabs(double x) + float sqrt(float x) + float fabs(float x) @cython.boundscheck(False) -cdef double sum_integral(np.ndarray[np.double_t, ndim=2, mode="c"] sat, +cdef float sum_integral(np.ndarray[float, ndim=2, mode="c"] sat, int r0, int c0, int r1, int c1): """ Using a summed area table / integral image, calculate the sum @@ -25,7 +25,7 @@ cdef double sum_integral(np.ndarray[np.double_t, ndim=2, mode="c"] sat, Parameters ---------- - sat : ndarray of double_t + sat : ndarray of float Summed area table / integral image. r0, c0 : int Top-left corner of block to be summed. @@ -37,7 +37,7 @@ cdef double sum_integral(np.ndarray[np.double_t, ndim=2, mode="c"] sat, S : int Sum over the given window. """ - cdef double S = 0 + cdef float S = 0 S += sat[r1, c1] @@ -53,27 +53,28 @@ cdef double sum_integral(np.ndarray[np.double_t, ndim=2, mode="c"] sat, @cython.boundscheck(False) -def match_template(np.ndarray[np.double_t, ndim=2, mode="c"] image, - np.ndarray[np.double_t, ndim=2, mode="c"] template, +def match_template(np.ndarray[float, ndim=2, mode="c"] image, + np.ndarray[float, ndim=2, mode="c"] template, int num_type): # convolve the image with template by frequency domain multiplication - cdef np.ndarray[np.double_t, ndim=2] result + 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.double) + mode="valid"), dtype=np.float32) # calculate squared integral images used for normalization - cdef np.ndarray[np.double_t, ndim=2, mode="c"] integral_sum - cdef np.ndarray[np.double_t, ndim=2, mode="c"] integral_sqr + cdef np.ndarray[float, ndim=2, mode="c"] integral_sum + cdef np.ndarray[float, ndim=2, mode="c"] integral_sqr if num_type == 1: integral_sum = integral.integral_image(image) integral_sqr = integral.integral_image(image**2) # use inversed area for accuracy - cdef double inv_area = 1.0 / (template.shape[0] * template.shape[1]) + 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 double template_norm - cdef double template_mean = np.mean(template) + cdef float template_norm + cdef float template_mean = np.mean(template) if num_type == 0: template_norm = sqrt((np.std(template) ** 2 + @@ -83,7 +84,7 @@ def match_template(np.ndarray[np.double_t, ndim=2, mode="c"] image, # define window of template size in squared integral image cdef int i, j - cdef double num, window_sum2, window_mean2, normed, t, + cdef float num, window_sum2, window_mean2, normed, t, # move window through convolution results, normalizing in the process for i in range(result.shape[0] - 1): for j in range(result.shape[1] - 1): diff --git a/skimage/feature/tests/test_template.py b/skimage/feature/tests/test_template.py index e92672da..33171abf 100644 --- a/skimage/feature/tests/test_template.py +++ b/skimage/feature/tests/test_template.py @@ -5,8 +5,9 @@ from numpy.random import randn def test_template(): size = 100 - image = np.zeros((400, 400)) + image = np.zeros((400, 400), dtype=np.float32) target = np.tri(size) + np.tri(size)[::-1] + target = target.astype(np.float32) target_positions = [(50, 50), (200, 200)] for x, y in target_positions: image[x:x + size, y:y + size] = target From 9d3072d3f8460d1d64a8ec1354461f03a45cae20 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 4 Feb 2012 00:34:34 -0500 Subject: [PATCH 20/48] Fix type errors in `match_template` by casting inputs to float32. --- skimage/feature/template.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/skimage/feature/template.py b/skimage/feature/template.py index 6bf02896..448326f1 100644 --- a/skimage/feature/template.py +++ b/skimage/feature/template.py @@ -3,6 +3,8 @@ import numpy as np import _template +from skimage.util.dtype import _convert + try: import cv opencv_available = True @@ -80,5 +82,7 @@ def match_template(image, template, method="norm-coeff"): method_num = 1 else: 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) From 425a4ea7074fa516a0355e44f5445a441eeabd2f Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 4 Feb 2012 01:35:27 -0500 Subject: [PATCH 21/48] Use `feature.peak_local_max` instead of custom peak detection. --- doc/examples/plot_template.py | 26 ++++++---------- skimage/feature/tests/test_template.py | 43 +++++++++++--------------- 2 files changed, 28 insertions(+), 41 deletions(-) diff --git a/doc/examples/plot_template.py b/doc/examples/plot_template.py index 3690b18e..aa8a1caa 100644 --- a/doc/examples/plot_template.py +++ b/doc/examples/plot_template.py @@ -13,7 +13,7 @@ a simple peak extraction algorithm. """ import numpy as np -from skimage.feature import match_template +from skimage.feature import match_template, peak_local_max from numpy.random import randn import matplotlib.pyplot as plt @@ -30,21 +30,14 @@ image += randn(400, 400)*2 # Match the template. result = match_template(image, target, method='norm-corr') -# peak extraction algorithm. -delta = 5 -found_positions = [] -for i in range(50): - index = np.argmax(result) - y, x = np.unravel_index(index, result.shape) - if not found_positions: - found_positions.append((x, y)) - for position in found_positions: - distance = np.sqrt((x - position[0]) ** 2 + (y - position[1]) ** 2) - if distance > delta: - found_positions.append((x, y)) - result[y, x] = 0 - if len(found_positions) == len(target_positions): - break +found_positions = peak_local_max(result) + +if len(found_positions) > 2: + # Keep the two maximum peaks. + intensities = result[tuple(found_positions.T)] + i_maxsort = np.argsort(intensities)[::-1] + found_positions = found_positions[i_maxsort][:2] + x_found, y_found = np.transpose(found_positions) plt.gray() @@ -67,3 +60,4 @@ plt.autoscale(tight=True) plt.axis('off') plt.show() + diff --git a/skimage/feature/tests/test_template.py b/skimage/feature/tests/test_template.py index 33171abf..0ec48f09 100644 --- a/skimage/feature/tests/test_template.py +++ b/skimage/feature/tests/test_template.py @@ -1,10 +1,13 @@ import numpy as np -from skimage.feature import match_template 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 def test_template(): size = 100 + # Type conversion of image and target not required but prevents warnings. image = np.zeros((400, 400), dtype=np.float32) target = np.tri(size) + np.tri(size)[::-1] target = target.astype(np.float32) @@ -16,31 +19,21 @@ def test_template(): for method in ["norm-corr", "norm-coeff"]: result = match_template(image, target, method=method) delta = 5 - found_positions = [] - # find the targets - for i in range(50): - index = np.argmax(result) - y, x = np.unravel_index(index, result.shape) - if not found_positions: - found_positions.append((x, y)) - for position in found_positions: - distance = np.sqrt((x - position[0]) ** 2 + - (y - position[1]) ** 2) - if distance > delta: - found_positions.append((x, y)) - result[y, x] = 0 - if len(found_positions) == len(target_positions): - break - for x, y in target_positions: - print x, y - found = False - for position in found_positions: - distance = np.sqrt((x - position[0]) ** 2 + - (y - position[1]) ** 2) - if distance < delta: - found = True - assert found + 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] + + # 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 + if __name__ == "__main__": from numpy import testing From 8c816109c960bc2125aa8499f121c13d68e2ef98 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 4 Feb 2012 09:10:43 -0500 Subject: [PATCH 22/48] Remove OpenCV version of match_template --- skimage/feature/template.py | 40 ------------------------------------- 1 file changed, 40 deletions(-) diff --git a/skimage/feature/template.py b/skimage/feature/template.py index 448326f1..b7ff70e1 100644 --- a/skimage/feature/template.py +++ b/skimage/feature/template.py @@ -5,46 +5,6 @@ import _template from skimage.util.dtype import _convert -try: - import cv - opencv_available = True -except ImportError: - opencv_available = False - - - -#XXX add to opencv backend once backend system in place -def match_template_cv(image, template, out=None, method="norm-coeff"): - """Finds a template in an image using normalized correlation. - - Parameters - ---------- - image : array_like, dtype=float - Image to process. - template : array_like, dtype=float - Template to locate. - out: array_like, dtype=float, optional - Optional destination. - Returns - ------- - output : ndarray, dtype=float - Correlation results between 0.0 and 1.0, maximum indicating the most - probable match. - """ - if not opencv_available: - raise ImportError("Opencv 2.0+ required") - if out == None: - out = np.empty((image.shape[0] - template.shape[0] + 1, - image.shape[1] - template.shape[1] + 1), - dtype=image.dtype) - if method == "norm-corr": - cv.MatchTemplate(image, template, out, cv.CV_TM_CCORR_NORMED) - elif method == "norm-coeff": - cv.MatchTemplate(image, template, out, cv.CV_TM_CCOEFF_NORMED) - else: - raise ValueError("Unknown template method: %s" % method) - return out - def match_template(image, template, method="norm-coeff"): """Finds a template in an image using normalized correlation. From 545bdab985a0c56f041b501601397e2586dfa694 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 4 Feb 2012 13:23:56 -0500 Subject: [PATCH 23/48] 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) From 3c7ae849b91c5320bc7c58cbadcc77c9a5a9e125 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 25 Feb 2012 23:54:43 -0500 Subject: [PATCH 24/48] Remove unnecessary truncation at boundary. --- skimage/feature/_template.pyx | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/skimage/feature/_template.pyx b/skimage/feature/_template.pyx index 44d78a47..77eea3f9 100644 --- a/skimage/feature/_template.pyx +++ b/skimage/feature/_template.pyx @@ -91,8 +91,8 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image, cdef int i, j cdef float num, window_sum2, window_mean2, normed, t, # move window through convolution results, normalizing in the process - for i in range(result.shape[0] - 1): - for j in range(result.shape[1] - 1): + for i in range(result.shape[0]): + for j in range(result.shape[1]): num = result[i, j] i_end = i + template.shape[0] j_end = j + template.shape[1] @@ -117,8 +117,5 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image, else: num = 0 result[i, j] = num - # zero boundaries - result[:, -1] = 0 - result[-1, :] = 0 return result From cbdea0d36e2eec31da2d761f70b73530e3dc2da2 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sun, 26 Feb 2012 00:20:44 -0500 Subject: [PATCH 25/48] Add `pad_output` argmument to `match_template`. --- skimage/feature/template.py | 31 +++++++++++++++++++++++-------- 1 file changed, 23 insertions(+), 8 deletions(-) diff --git a/skimage/feature/template.py b/skimage/feature/template.py index 40306fab..643e67c6 100644 --- a/skimage/feature/template.py +++ b/skimage/feature/template.py @@ -6,7 +6,7 @@ import _template from skimage.util.dtype import _convert -def match_template(image, template, method='norm-coeff'): +def match_template(image, template, method='norm-coeff', pad_output=True): """Finds a template in an image using normalized correlation. TODO: The output is currently smaller than the input image due to @@ -14,11 +14,11 @@ def match_template(image, template, method='norm-coeff'): Parameters ---------- - image : array_like, dtype=float + image : array_like Image to process. - template : array_like, dtype=float + template : array_like Template to locate. - method: str (default 'norm-coeff') + method : str The correlation method used in scanning. 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. @@ -32,17 +32,32 @@ def match_template(image, template, method='norm-coeff'): 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) * Sum(X',Y')[I(x + X', y + Y')] + 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)` + for an `(M, N)` image and an `(m, n)` template. Returns ------- - output : ndarray, dtype=float - Correlation results between 0.0 and 1.0, maximum indicating the most - probable match. + 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. """ 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) + result = _template.match_template(image, template, method) + + if pad_output: + h, w = result.shape + full_result = np.zeros(image.shape, dtype=np.float32) + full_result[:h, :w] = result + return full_result + else: + return result From a87bcb2d7325925ddf2c68da0aab37d3cd652c55 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sun, 26 Feb 2012 00:21:47 -0500 Subject: [PATCH 26/48] DOC: demonstrate where the template is matched in an image. --- doc/examples/plot_template.py | 33 ++++++++++++++++++--------------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/doc/examples/plot_template.py b/doc/examples/plot_template.py index aa8a1caa..fe6549f6 100644 --- a/doc/examples/plot_template.py +++ b/doc/examples/plot_template.py @@ -9,7 +9,9 @@ techniques to find instances of the "target image" in the "test image". The output of ``match_template`` is an image where we can easily identify peaks by eye. We mark the locations of matches (red dots), which are detected using -a simple peak extraction algorithm. +a simple peak extraction algorithm. Note that the peaks in the output of +``match_template`` correspond to the origin (i.e. top-left corner) of the +template. """ import numpy as np @@ -40,24 +42,25 @@ if len(found_positions) > 2: x_found, y_found = np.transpose(found_positions) + +fig, (ax0, ax1, ax2) = plt.subplots(ncols=3, figsize=(8, 3)) plt.gray() -plt.subplot(1, 3, 1) -plt.imshow(target) -plt.title("Target image") -plt.axis('off') +ax0.imshow(target) +ax0.set_title("Target image") -plt.subplot(1, 3, 2) -plt.imshow(image) -plt.title("Test image") -plt.axis('off') +ax1.imshow(image) +ax1.plot(x_found, y_found, 'ro', alpha=0.5) +ax1.set_title("Test image") +ax1.autoscale(tight=True) -plt.subplot(1, 3, 3) -plt.imshow(result) -plt.plot(x_found, y_found, 'ro') -plt.title("Result from\n``match_template``") -plt.autoscale(tight=True) -plt.axis('off') +ax2.imshow(result) +ax2.plot(x_found, y_found, 'ro', alpha=0.5) +ax2.set_title("Result from\n``match_template``") +ax2.autoscale(tight=True) + +for ax in (ax0, ax1, ax2): + ax.axis('off') plt.show() From 2e8fcef89b1cdd5c0910a0733f33ffb6c339f7a3 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sun, 26 Feb 2012 10:22:54 -0500 Subject: [PATCH 27/48] Simply equation in docstring of `match_docstring`. --- skimage/feature/template.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/skimage/feature/template.py b/skimage/feature/template.py index 643e67c6..39713e43 100644 --- a/skimage/feature/template.py +++ b/skimage/feature/template.py @@ -21,17 +21,20 @@ def match_template(image, template, method='norm-coeff', pad_output=True): method : str The correlation method used in scanning. 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. + All sums are done over X = 0..w-1 and Y = 0..h-1 of the template. 'norm-coeff': - 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]) + 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(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]) + 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) - 1/(w * h) * Sum(X',Y')[T(X', Y')] - I'(x + X, y + Y) = I(x + X, y + Y) - - 1/(w * h) * Sum(X',Y')[I(x + X', y + Y')] + + 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)` From 883dd24cdb4b518d6b5bd30813afe4d42541674d Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Wed, 29 Feb 2012 00:34:04 -0500 Subject: [PATCH 28/48] Fix bug when indexing into summed-area table. --- skimage/feature/_template.pyx | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/skimage/feature/_template.pyx b/skimage/feature/_template.pyx index 77eea3f9..0c3145f1 100644 --- a/skimage/feature/_template.pyx +++ b/skimage/feature/_template.pyx @@ -94,8 +94,10 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image, for i in range(result.shape[0]): for j in range(result.shape[1]): num = result[i, j] - i_end = i + template.shape[0] - j_end = j + template.shape[1] + # subtract 1 because `i_end` and `j_end` are used for indexing into + # summed-area table, instead of slicing windows of the image. + i_end = i + template.shape[0] - 1 + j_end = j + template.shape[1] - 1 window_mean2 = 0 if method == 'norm-coeff': From 753e999a7ab76561d1ec7ddee9124dac4b05174b Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Wed, 29 Feb 2012 00:54:37 -0500 Subject: [PATCH 29/48] Minor cleanup * Rename some variables. * Delete some old comments. * Move some variable initializations to the top of the function. --- skimage/feature/_template.pyx | 39 ++++++++++++++++------------------- 1 file changed, 18 insertions(+), 21 deletions(-) diff --git a/skimage/feature/_template.pyx b/skimage/feature/_template.pyx index 0c3145f1..9052fd3c 100644 --- a/skimage/feature/_template.pyx +++ b/skimage/feature/_template.pyx @@ -56,24 +56,22 @@ cdef float sum_integral(np.ndarray[float, ndim=2, mode="c"] sat, def match_template(np.ndarray[float, ndim=2, mode="c"] image, np.ndarray[float, ndim=2, mode="c"] template, str method): - # convolve the image with template by frequency domain multiplication 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 + cdef float template_mean = np.mean(template) + cdef float template_norm + cdef float inv_area + # 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 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]) - cdef float template_norm - cdef float template_mean = np.mean(template) + inv_area = 1.0 / (template.shape[0] * template.shape[1]) if method == 'norm-corr': # calculate template norm according to the following: @@ -87,9 +85,8 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image, else: template_norm = sqrt((template_mean ** 2)) / sqrt(inv_area) - # define window of template size in squared integral image cdef int i, j - cdef float num, window_sum2, window_mean2, normed, t, + cdef float num, den, window_sqr_sum, window_mean_sqr, window_sum, # move window through convolution results, normalizing in the process for i in range(result.shape[0]): for j in range(result.shape[1]): @@ -99,19 +96,19 @@ 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_mean2 = 0 + window_mean_sqr = 0 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_end, j_end) + 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 - normed = sqrt(window_sum2 - window_mean2) * template_norm + window_sqr_sum = sum_integral(integral_sqr, i, j, i_end, j_end) + + den = sqrt(window_sqr_sum - window_mean_sqr) * template_norm # enforce some limits - if fabs(num) < normed: - num /= normed - elif fabs(num) < normed*1.125: + if fabs(num) < den: + num /= den + elif fabs(num) < den * 1.125: if num > 0: num = 1 else: From 5682d27eb089108bb3c106e22ae6a7ce6baa0c91 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Thu, 1 Mar 2012 23:00:34 -0500 Subject: [PATCH 30/48] 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__": From 94262bc8d6d12ae6910e7e1af33b963ccc296b61 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Thu, 1 Mar 2012 23:05:17 -0500 Subject: [PATCH 31/48] Tweak example to reduce confusion. Set background 1 so that the template has values equally-above and -below the mean image value. This change makes the template results left-right symmetric, which may reduce confusion. --- doc/examples/plot_template.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/examples/plot_template.py b/doc/examples/plot_template.py index 3aea93d9..3ebb9070 100644 --- a/doc/examples/plot_template.py +++ b/doc/examples/plot_template.py @@ -16,18 +16,18 @@ template. import numpy as np from skimage.feature import match_template, peak_local_max -from numpy.random import randn import matplotlib.pyplot as plt # We first construct a simple image target: size = 100 target = np.tri(size) + np.tri(size)[::-1] # place target in an image at two positions, and add noise. -image = np.zeros((400, 400)) +image = np.ones((400, 400)) 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 result = match_template(image, target) From f21f032bfee339bf8e1b422a0627b233dc3e7508 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Thu, 1 Mar 2012 23:23:40 -0500 Subject: [PATCH 32/48] Add normalization test and remove unnecessary clipping. --- skimage/feature/_template.pyx | 12 +++------ skimage/feature/tests/test_template.py | 36 ++++++++++++++++++++++++++ 2 files changed, 39 insertions(+), 9 deletions(-) diff --git a/skimage/feature/_template.pyx b/skimage/feature/_template.pyx index 63ffa3d9..9e9a2148 100644 --- a/skimage/feature/_template.pyx +++ b/skimage/feature/_template.pyx @@ -120,16 +120,10 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image, window_sqr_sum = sum_integral(integral_sqr, i, j, i_end, j_end) den = sqrt((window_sqr_sum - window_mean_sqr) * template_ssd) - # enforce some limits - if fabs(num) < den: - num /= den - elif fabs(num) < den * 1.125: - if num > 0: - num = 1 - else: - num = -1 - else: + if den == 0: num = 0 + else: + num /= den result[i, j] = num return result diff --git a/skimage/feature/tests/test_template.py b/skimage/feature/tests/test_template.py index a3af9954..36b72fd5 100644 --- a/skimage/feature/tests/test_template.py +++ b/skimage/feature/tests/test_template.py @@ -34,6 +34,42 @@ def test_template(): yield assert_close, xy, xy_target +def test_normalization(): + """Test that `match_template` gives the correct normalization. + + Normalization gives 1 for a perfect match and -1 for an inverted-match. + This test adds positive and negative squares to a zero-array and matches + the array with a positive template. + """ + n = 5 + N = 20 + ipos, jpos = (2, 3) + ineg, jneg = (12, 11) + image = np.zeros((N, N)) + image[ipos:ipos + n, jpos:jpos + n] = 10 + image[ineg:ineg + n, jneg:jneg + n] = -10 + + # white square with a black border + template = np.zeros((n+2, n+2)) + template[1:1+n, 1:1+n] = 1 + + result = match_template(image, template) + + # get the max and min results. + sorted_result = np.argsort(result.flat) + iflat_min = sorted_result[0] + iflat_max = sorted_result[-1] + min_result = np.unravel_index(iflat_min, (N, N)) + max_result = np.unravel_index(iflat_max, (N, N)) + + # shift result by 1 because of template border + assert np.all((np.array(min_result) + 1) == (ineg, jneg)) + assert np.all((np.array(max_result) + 1) == (ipos, jpos)) + + assert np.allclose(result.flat[iflat_min], -1) + assert np.allclose(result.flat[iflat_max], 1) + + if __name__ == "__main__": from numpy import testing testing.run_module_suite() From 7caf85a5c1efe9ae17f6aa62c3636b489a33ffdb Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Thu, 1 Mar 2012 23:32:48 -0500 Subject: [PATCH 33/48] Renaming for clarity. --- skimage/feature/_template.pyx | 34 ++++++++++++++++------------------ 1 file changed, 16 insertions(+), 18 deletions(-) diff --git a/skimage/feature/_template.pyx b/skimage/feature/_template.pyx index 9e9a2148..b0f6dc66 100644 --- a/skimage/feature/_template.pyx +++ b/skimage/feature/_template.pyx @@ -43,7 +43,7 @@ cdef extern from "math.h": @cython.boundscheck(False) -cdef float sum_integral(np.ndarray[float, ndim=2, mode="c"] sat, +cdef float integrate(np.ndarray[float, ndim=2, mode="c"] sat, int r0, int c0, int r1, int c1): """ Using a summed area table / integral image, calculate the sum @@ -85,15 +85,15 @@ 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): - 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 np.ndarray[float, ndim=2, mode="c"] corr + cdef np.ndarray[float, ndim=2, mode="c"] image_sat + cdef np.ndarray[float, ndim=2, mode="c"] image_sqr_sat cdef float template_mean = np.mean(template) cdef float template_ssd cdef float inv_area - integral_sum = integral.integral_image(image) - integral_sqr = integral.integral_image(image**2) + image_sat = integral.integral_image(image) + image_sqr_sat = integral.integral_image(image**2) template -= template_mean template_ssd = np.sum(template**2) @@ -101,29 +101,27 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image, 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) + corr = np.ascontiguousarray(fftconvolve(image, np.fliplr(template), + mode="valid"), dtype=np.float32) cdef int i, j - cdef float num, den, window_sqr_sum, window_mean_sqr, window_sum, + cdef float den, window_sqr_sum, window_mean_sqr, window_sum, # move window through convolution results, normalizing in the process - for i in range(result.shape[0]): - for j in range(result.shape[1]): - num = result[i, j] + for i in range(corr.shape[0]): + for j in range(corr.shape[1]): # subtract 1 because `i_end` and `j_end` are used for indexing into # summed-area table, instead of slicing windows of the image. i_end = i + template.shape[0] - 1 j_end = j + template.shape[1] - 1 - window_sum = sum_integral(integral_sum, i, j, i_end, j_end) + window_sum = integrate(image_sat, 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) + window_sqr_sum = integrate(image_sqr_sat, i, j, i_end, j_end) den = sqrt((window_sqr_sum - window_mean_sqr) * template_ssd) if den == 0: - num = 0 + corr[i, j] = 0 else: - num /= den - result[i, j] = num - return result + corr[i, j] /= den + return corr From 5d79537566c0e1df94326e23ad5bc4162ce27694 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 24 Mar 2012 17:42:36 -0400 Subject: [PATCH 34/48] Fix convolution input to get correlation results. --- skimage/feature/_template.pyx | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/skimage/feature/_template.pyx b/skimage/feature/_template.pyx index b0f6dc66..7f31d31e 100644 --- a/skimage/feature/_template.pyx +++ b/skimage/feature/_template.pyx @@ -101,8 +101,10 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image, inv_area = 1.0 / (template.shape[0] * template.shape[1]) # when `dtype=float` is used, ascontiguousarray returns ``double``. - corr = np.ascontiguousarray(fftconvolve(image, np.fliplr(template), - mode="valid"), dtype=np.float32) + corr = np.ascontiguousarray(fftconvolve(image, + np.flipud(np.fliplr(template)), + mode="valid"), + dtype=np.float32) cdef int i, j cdef float den, window_sqr_sum, window_mean_sqr, window_sum, From 43ed4b9fb762a99637a714a11837284c85d8f74b Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Sat, 24 Mar 2012 20:06:24 -0400 Subject: [PATCH 35/48] Pad input image instead of output. --- skimage/feature/template.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/skimage/feature/template.py b/skimage/feature/template.py index 3c85985f..398e7c78 100644 --- a/skimage/feature/template.py +++ b/skimage/feature/template.py @@ -35,13 +35,16 @@ def match_template(image, template, pad_output=True): """ image = _convert(image, np.float32) template = _convert(template, np.float32) - result = _template.match_template(image, template) if pad_output: - h, w = result.shape - full_result = np.zeros(image.shape, dtype=np.float32) - full_result[:h, :w] = result - return full_result - else: - return result + 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 From b2006da3461106ba41834acb6d3c81e45f4ea08e Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Mon, 26 Mar 2012 20:17:55 -0400 Subject: [PATCH 36/48] Change `flipud`/`fliplr` to array indexing. --- skimage/feature/_template.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/feature/_template.pyx b/skimage/feature/_template.pyx index 7f31d31e..4fa776d4 100644 --- a/skimage/feature/_template.pyx +++ b/skimage/feature/_template.pyx @@ -102,7 +102,7 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image, # when `dtype=float` is used, ascontiguousarray returns ``double``. corr = np.ascontiguousarray(fftconvolve(image, - np.flipud(np.fliplr(template)), + template[::-1, ::-1], mode="valid"), dtype=np.float32) From 5e49eb4fb6bce9cdeae515590530b78e4dde89d9 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Mon, 26 Mar 2012 20:26:44 -0400 Subject: [PATCH 37/48] Add alternate example for `match_template`. --- doc/examples/plot_match_face_template.py | 41 ++++++++++++++++++++++++ 1 file changed, 41 insertions(+) create mode 100644 doc/examples/plot_match_face_template.py diff --git a/doc/examples/plot_match_face_template.py b/doc/examples/plot_match_face_template.py new file mode 100644 index 00000000..898870e5 --- /dev/null +++ b/doc/examples/plot_match_face_template.py @@ -0,0 +1,41 @@ +""" +================= +Template Matching +================= + +In this example, we use template matching to identify the occurrence of an +image patch (in this case, a sub-image centered on the camera man's head). +Since there's only a single match, the maximum value in the `match_template` +result` corresponds to the head location. If you expect multiple matches, you +should use a proper peak-finding function. + +""" + +import numpy as np +import matplotlib.pyplot as plt +from skimage import data +from skimage.feature import match_template + +image = data.camera() +head = image[70:170, 180:280] + +result = match_template(image, head) + +fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(8, 4)) + +ax1.imshow(head) +ax1.set_axis_off() +ax1.set_title('template') + +ax2.imshow(image) +ax2.set_axis_off() +ax2.set_title('image') + +# highlight matched region +xy = np.unravel_index(np.argmax(result), image.shape)[::-1] # -1 flips ij to xy +wface, hface = head.shape +rect = plt.Rectangle(xy, wface, hface, edgecolor='r', facecolor='none') +ax2.add_patch(rect) + +plt.show() + From 0c1e3541b37e0c46b7de6747709f915f7fa1a727 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Mon, 26 Mar 2012 20:43:45 -0400 Subject: [PATCH 38/48] Check that image is larger than template. --- skimage/feature/template.py | 2 ++ skimage/feature/tests/test_template.py | 6 ++++++ 2 files changed, 8 insertions(+) diff --git a/skimage/feature/template.py b/skimage/feature/template.py index 398e7c78..67064d5a 100644 --- a/skimage/feature/template.py +++ b/skimage/feature/template.py @@ -33,6 +33,8 @@ def match_template(image, template, pad_output=True): beyond the image edges. """ + if np.any(np.less(image.shape, template.shape)): + raise ValueError("Image must be larger than template.") image = _convert(image, np.float32) template = _convert(template, np.float32) diff --git a/skimage/feature/tests/test_template.py b/skimage/feature/tests/test_template.py index 36b72fd5..5ad1fb9b 100644 --- a/skimage/feature/tests/test_template.py +++ b/skimage/feature/tests/test_template.py @@ -70,6 +70,12 @@ def test_normalization(): assert np.allclose(result.flat[iflat_max], 1) +def test_switched_arguments(): + image = np.ones((5, 5)) + template = np.ones((3, 3)) + np.testing.assert_raises(ValueError, match_template, template, image) + + if __name__ == "__main__": from numpy import testing testing.run_module_suite() From ad99285bc070519644c2485ef75e048ac6614c4f Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Mon, 26 Mar 2012 22:27:42 -0400 Subject: [PATCH 39/48] Prevent `match_template` from returning NaNs. --- skimage/feature/_template.pyx | 10 +++++----- skimage/feature/tests/test_template.py | 15 +++++++++++++++ 2 files changed, 20 insertions(+), 5 deletions(-) diff --git a/skimage/feature/_template.pyx b/skimage/feature/_template.pyx index 4fa776d4..b83761a8 100644 --- a/skimage/feature/_template.pyx +++ b/skimage/feature/_template.pyx @@ -119,11 +119,11 @@ def match_template(np.ndarray[float, ndim=2, mode="c"] image, window_sum = integrate(image_sat, i, j, i_end, j_end) window_mean_sqr = window_sum * window_sum * inv_area window_sqr_sum = integrate(image_sqr_sat, i, j, i_end, j_end) - den = sqrt((window_sqr_sum - window_mean_sqr) * template_ssd) - - if den == 0: + if window_sqr_sum <= window_mean_sqr: corr[i, j] = 0 - else: - corr[i, j] /= den + continue + + den = sqrt((window_sqr_sum - window_mean_sqr) * template_ssd) + corr[i, j] /= den return corr diff --git a/skimage/feature/tests/test_template.py b/skimage/feature/tests/test_template.py index 5ad1fb9b..dd24cbd6 100644 --- a/skimage/feature/tests/test_template.py +++ b/skimage/feature/tests/test_template.py @@ -70,6 +70,21 @@ def test_normalization(): assert np.allclose(result.flat[iflat_max], 1) +def test_no_nans(): + """Test that `match_template` doesn't return NaN values. + + When image values are only slightly different, floating-point errors can + cause a subtraction inside of a square root to go negative (without an + explicit check that was added to `match_template`). + """ + np.random.seed(1) + image = 10000 + np.random.normal(size=(20, 20)) + template = np.ones((6, 6)) + template[:3, :] = 0 + result = match_template(image, template) + assert not np.any(np.isnan(result)) + + def test_switched_arguments(): image = np.ones((5, 5)) template = np.ones((3, 3)) From 383ca6220a4667230c652392b840fa82b18eea7a Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Thu, 29 Mar 2012 21:49:27 -0400 Subject: [PATCH 40/48] Set `match_template` to default to no padding and fix test. --- skimage/feature/template.py | 2 +- skimage/feature/tests/test_template.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/skimage/feature/template.py b/skimage/feature/template.py index 67064d5a..b328131f 100644 --- a/skimage/feature/template.py +++ b/skimage/feature/template.py @@ -6,7 +6,7 @@ import _template from skimage.util.dtype import _convert -def match_template(image, template, pad_output=True): +def match_template(image, template, pad_output=False): """Match a template to an image using normalized correlation. The output is an array with values between -1.0 and 1.0, which correspond diff --git a/skimage/feature/tests/test_template.py b/skimage/feature/tests/test_template.py index dd24cbd6..c6df75d0 100644 --- a/skimage/feature/tests/test_template.py +++ b/skimage/feature/tests/test_template.py @@ -59,8 +59,8 @@ def test_normalization(): sorted_result = np.argsort(result.flat) iflat_min = sorted_result[0] iflat_max = sorted_result[-1] - min_result = np.unravel_index(iflat_min, (N, N)) - max_result = np.unravel_index(iflat_max, (N, N)) + min_result = np.unravel_index(iflat_min, result.shape) + max_result = np.unravel_index(iflat_max, result.shape) # shift result by 1 because of template border assert np.all((np.array(min_result) + 1) == (ineg, jneg)) From 193d5abd3c49d25da61cf4385f39f967e7c2a598 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Thu, 29 Mar 2012 23:58:33 -0400 Subject: [PATCH 41/48] Rename `pad_output` parameter to `pad_input`. --- skimage/feature/template.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/skimage/feature/template.py b/skimage/feature/template.py index b328131f..6ad45926 100644 --- a/skimage/feature/template.py +++ b/skimage/feature/template.py @@ -6,12 +6,11 @@ import _template from skimage.util.dtype import _convert -def match_template(image, template, pad_output=False): +def match_template(image, template, pad_input=False): """Match a template to an 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's *origin* (i.e. its top-left - corner) is found at that position. + to the probability that the template is found at that position. Parameters ---------- @@ -19,18 +18,19 @@ def match_template(image, template, pad_output=False): Image to process. template : array_like Template to locate. - pad_output : bool - If True, pad output array to be the same size as the input image. + pad_input : bool + If True, pad `image` with image mean so that output is the same size as + the image, and output values correspond to the template center. Otherwise, the output is an array with shape `(M - m + 1, N - n + 1)` - for an `(M, N)` image and an `(m, n)` template. + for an `(M, N)` image and an `(m, n)` template, and matches correspond + to origin (top-left corner) of the template. Returns ------- output : ndarray - 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. + Correlation results between -1.0 and 1.0. For an `(M, N)` image and an + `(m, n)` template, the `output` is `(M - m + 1, N - n + 1)` when + `pad_input = False` and `(M, N)` when `pad_input = True`. """ if np.any(np.less(image.shape, template.shape)): @@ -38,7 +38,7 @@ def match_template(image, template, pad_output=False): image = _convert(image, np.float32) template = _convert(template, np.float32) - if pad_output: + 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 From a2f91585ef5772d5c08eb50bcbb1841f0aabf494 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Thu, 29 Mar 2012 23:59:09 -0400 Subject: [PATCH 42/48] Add test for `pad_input = True`. --- skimage/feature/tests/test_template.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/skimage/feature/tests/test_template.py b/skimage/feature/tests/test_template.py index c6df75d0..2eedb922 100644 --- a/skimage/feature/tests/test_template.py +++ b/skimage/feature/tests/test_template.py @@ -1,6 +1,7 @@ import numpy as np from numpy.testing import assert_array_almost_equal as assert_close +from skimage.morphology import diamond from skimage.feature import match_template, peak_local_max @@ -91,6 +92,26 @@ def test_switched_arguments(): np.testing.assert_raises(ValueError, match_template, template, image) +def test_pad_input(): + template = 10.0 * diamond(2) + + image = np.zeros((9, 19)) + mid = slice(2, 7) + image[mid, :3] = -template[:, -3:] # half min template centered at 0 + image[mid, 4:9] = template # full max template centered at 6 + 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) + + # get the max and min results. + sorted_result = np.argsort(result.flat) + i, j = np.unravel_index(sorted_result[:2], result.shape) + assert_close(j, (12, 0)) + i, j = np.unravel_index(sorted_result[-2:], result.shape) + assert_close(j, (18, 6)) + + if __name__ == "__main__": from numpy import testing testing.run_module_suite() From 06bd3ebd6fc8807be35ddd41dc88cc9a42c6faa4 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Fri, 30 Mar 2012 00:10:03 -0400 Subject: [PATCH 43/48] Add doctest example to `match_template`. --- skimage/feature/template.py | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/skimage/feature/template.py b/skimage/feature/template.py index 6ad45926..e5ba8c21 100644 --- a/skimage/feature/template.py +++ b/skimage/feature/template.py @@ -32,6 +32,38 @@ 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`. + Examples + -------- + >>> template = np.zeros((3, 3)) + >>> template[1, 1] = 1 + >>> print template + [[ 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.]] + >>> 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. ]] + >>> 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]] """ if np.any(np.less(image.shape, template.shape)): raise ValueError("Image must be larger than template.") From 4c3f8a36bad746218a16e29862c0466121be4646 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Fri, 30 Mar 2012 00:32:33 -0400 Subject: [PATCH 44/48] Fix template example. --- doc/examples/plot_match_face_template.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/examples/plot_match_face_template.py b/doc/examples/plot_match_face_template.py index 898870e5..49905941 100644 --- a/doc/examples/plot_match_face_template.py +++ b/doc/examples/plot_match_face_template.py @@ -32,7 +32,7 @@ ax2.set_axis_off() ax2.set_title('image') # highlight matched region -xy = np.unravel_index(np.argmax(result), image.shape)[::-1] # -1 flips ij to xy +xy = np.unravel_index(np.argmax(result), result.shape)[::-1] #-1 flips ij to xy wface, hface = head.shape rect = plt.Rectangle(xy, wface, hface, edgecolor='r', facecolor='none') ax2.add_patch(rect) From 8dbf6f4c581430ae7393d1ed0c5f0b377ffebd7e Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Fri, 30 Mar 2012 10:18:03 -0400 Subject: [PATCH 45/48] Fix shape unpacking ((height, width), not (w, h)). --- doc/examples/plot_match_face_template.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/examples/plot_match_face_template.py b/doc/examples/plot_match_face_template.py index 49905941..d3cb7905 100644 --- a/doc/examples/plot_match_face_template.py +++ b/doc/examples/plot_match_face_template.py @@ -33,7 +33,7 @@ ax2.set_title('image') # highlight matched region xy = np.unravel_index(np.argmax(result), result.shape)[::-1] #-1 flips ij to xy -wface, hface = head.shape +hface, wface = head.shape rect = plt.Rectangle(xy, wface, hface, edgecolor='r', facecolor='none') ax2.add_patch(rect) From f3e91020f0426fedfe229e94bf1ddc69dd64a136 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Fri, 30 Mar 2012 10:24:16 -0400 Subject: [PATCH 46/48] Add new example plot for `match_template`. --- doc/examples/plot_template_alt.py | 56 +++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) create mode 100644 doc/examples/plot_template_alt.py diff --git a/doc/examples/plot_template_alt.py b/doc/examples/plot_template_alt.py new file mode 100644 index 00000000..65b4571b --- /dev/null +++ b/doc/examples/plot_template_alt.py @@ -0,0 +1,56 @@ +""" +================= +Template Matching +================= + +In this example, we use template matching to identify the occurrence of an +image patch (in this case, a sub-image centered on a single coin). Here, we +return a single match (the exact same coin), so the maximum value in the +``match_template`` result corresponds to the coin location. The other coins +look similar, and thus have local maxima; if you expect multiple matches, you +should use a proper peak-finding function. + +The ``match_template`` function uses fast, normalized cross-correlation [1]_ +to find instances of the template in the image. Note that the peaks in the +output of ``match_template`` correspond to the origin (i.e. top-left corner) of +the template. + +.. [1] J. P. Lewis, "Fast Normalized Cross-Correlation", Industrial Light and + Magic. +""" + +import numpy as np +import matplotlib.pyplot as plt +from skimage import data +from skimage.feature import match_template + +image = data.coins() +coin = image[170:220, 75:130] + +result = match_template(image, coin) +ij = np.unravel_index(np.argmax(result), result.shape) +x, y = ij[::-1] + +fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(8, 3)) + +ax1.imshow(coin) +ax1.set_axis_off() +ax1.set_title('template') + +ax2.imshow(image) +ax2.set_axis_off() +ax2.set_title('image') +# highlight matched region +hcoin, wcoin = coin.shape +rect = plt.Rectangle((x, y), wcoin, hcoin, edgecolor='r', facecolor='none') +ax2.add_patch(rect) + +ax3.imshow(result) +ax3.set_axis_off() +ax3.set_title('`match_template`\nresult') +# highlight matched region +ax3.autoscale(False) +ax3.plot(x, y, 'o', markeredgecolor='r', markerfacecolor='none', markersize=10) + +plt.show() + From 3c3c95b40631c80c746782c28177f5a2fde8d6f0 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Mon, 16 Apr 2012 20:16:08 -0400 Subject: [PATCH 47/48] DOC: Replace template example with alternate example. And remove other alternate example. --- doc/examples/plot_match_face_template.py | 41 ------------ doc/examples/plot_template.py | 81 +++++++++++------------- doc/examples/plot_template_alt.py | 56 ---------------- 3 files changed, 36 insertions(+), 142 deletions(-) delete mode 100644 doc/examples/plot_match_face_template.py delete mode 100644 doc/examples/plot_template_alt.py diff --git a/doc/examples/plot_match_face_template.py b/doc/examples/plot_match_face_template.py deleted file mode 100644 index d3cb7905..00000000 --- a/doc/examples/plot_match_face_template.py +++ /dev/null @@ -1,41 +0,0 @@ -""" -================= -Template Matching -================= - -In this example, we use template matching to identify the occurrence of an -image patch (in this case, a sub-image centered on the camera man's head). -Since there's only a single match, the maximum value in the `match_template` -result` corresponds to the head location. If you expect multiple matches, you -should use a proper peak-finding function. - -""" - -import numpy as np -import matplotlib.pyplot as plt -from skimage import data -from skimage.feature import match_template - -image = data.camera() -head = image[70:170, 180:280] - -result = match_template(image, head) - -fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(8, 4)) - -ax1.imshow(head) -ax1.set_axis_off() -ax1.set_title('template') - -ax2.imshow(image) -ax2.set_axis_off() -ax2.set_title('image') - -# highlight matched region -xy = np.unravel_index(np.argmax(result), result.shape)[::-1] #-1 flips ij to xy -hface, wface = head.shape -rect = plt.Rectangle(xy, wface, hface, edgecolor='r', facecolor='none') -ax2.add_patch(rect) - -plt.show() - diff --git a/doc/examples/plot_template.py b/doc/examples/plot_template.py index 3ebb9070..65b4571b 100644 --- a/doc/examples/plot_template.py +++ b/doc/examples/plot_template.py @@ -4,62 +4,53 @@ Template Matching ================= In this example, we use template matching to identify the occurrence of an -object in an image. The ``match_template`` function uses normalised correlation -techniques to find instances of the "target image" in the "test image". +image patch (in this case, a sub-image centered on a single coin). Here, we +return a single match (the exact same coin), so the maximum value in the +``match_template`` result corresponds to the coin location. The other coins +look similar, and thus have local maxima; if you expect multiple matches, you +should use a proper peak-finding function. -The output of ``match_template`` is an image where we can easily identify peaks -by eye. We mark the locations of matches (red dots), which are detected using -a simple peak extraction algorithm. Note that the peaks in the output of -``match_template`` correspond to the origin (i.e. top-left corner) of the -template. +The ``match_template`` function uses fast, normalized cross-correlation [1]_ +to find instances of the template in the image. Note that the peaks in the +output of ``match_template`` correspond to the origin (i.e. top-left corner) of +the template. + +.. [1] J. P. Lewis, "Fast Normalized Cross-Correlation", Industrial Light and + Magic. """ import numpy as np -from skimage.feature import match_template, peak_local_max import matplotlib.pyplot as plt +from skimage import data +from skimage.feature import match_template -# We first construct a simple image target: -size = 100 -target = np.tri(size) + np.tri(size)[::-1] -# place target in an image at two positions, and add noise. -image = np.ones((400, 400)) -target_positions = [(50, 50), (200, 200)] -for x, y in target_positions: - image[x:x+size, y:y+size] = target -np.random.seed(1) -image += np.random.randn(400, 400)*2 +image = data.coins() +coin = image[170:220, 75:130] -result = match_template(image, target) +result = match_template(image, coin) +ij = np.unravel_index(np.argmax(result), result.shape) +x, y = ij[::-1] -found_positions = peak_local_max(result) +fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(8, 3)) -if len(found_positions) > 2: - # Keep the two maximum peaks. - intensities = result[tuple(found_positions.T)] - i_maxsort = np.argsort(intensities)[::-1] - found_positions = found_positions[i_maxsort][:2] +ax1.imshow(coin) +ax1.set_axis_off() +ax1.set_title('template') -x_found, y_found = np.transpose(found_positions) +ax2.imshow(image) +ax2.set_axis_off() +ax2.set_title('image') +# highlight matched region +hcoin, wcoin = coin.shape +rect = plt.Rectangle((x, y), wcoin, hcoin, edgecolor='r', facecolor='none') +ax2.add_patch(rect) - -fig, (ax0, ax1, ax2) = plt.subplots(ncols=3, figsize=(8, 3)) -plt.gray() - -ax0.imshow(target) -ax0.set_title("Target image") - -ax1.imshow(image) -ax1.plot(x_found, y_found, 'ro', alpha=0.5) -ax1.set_title("Test image") -ax1.autoscale(tight=True) - -ax2.imshow(result) -ax2.plot(x_found, y_found, 'ro', alpha=0.5) -ax2.set_title("Result from\n``match_template``") -ax2.autoscale(tight=True) - -for ax in (ax0, ax1, ax2): - ax.axis('off') +ax3.imshow(result) +ax3.set_axis_off() +ax3.set_title('`match_template`\nresult') +# highlight matched region +ax3.autoscale(False) +ax3.plot(x, y, 'o', markeredgecolor='r', markerfacecolor='none', markersize=10) plt.show() diff --git a/doc/examples/plot_template_alt.py b/doc/examples/plot_template_alt.py deleted file mode 100644 index 65b4571b..00000000 --- a/doc/examples/plot_template_alt.py +++ /dev/null @@ -1,56 +0,0 @@ -""" -================= -Template Matching -================= - -In this example, we use template matching to identify the occurrence of an -image patch (in this case, a sub-image centered on a single coin). Here, we -return a single match (the exact same coin), so the maximum value in the -``match_template`` result corresponds to the coin location. The other coins -look similar, and thus have local maxima; if you expect multiple matches, you -should use a proper peak-finding function. - -The ``match_template`` function uses fast, normalized cross-correlation [1]_ -to find instances of the template in the image. Note that the peaks in the -output of ``match_template`` correspond to the origin (i.e. top-left corner) of -the template. - -.. [1] J. P. Lewis, "Fast Normalized Cross-Correlation", Industrial Light and - Magic. -""" - -import numpy as np -import matplotlib.pyplot as plt -from skimage import data -from skimage.feature import match_template - -image = data.coins() -coin = image[170:220, 75:130] - -result = match_template(image, coin) -ij = np.unravel_index(np.argmax(result), result.shape) -x, y = ij[::-1] - -fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, figsize=(8, 3)) - -ax1.imshow(coin) -ax1.set_axis_off() -ax1.set_title('template') - -ax2.imshow(image) -ax2.set_axis_off() -ax2.set_title('image') -# highlight matched region -hcoin, wcoin = coin.shape -rect = plt.Rectangle((x, y), wcoin, hcoin, edgecolor='r', facecolor='none') -ax2.add_patch(rect) - -ax3.imshow(result) -ax3.set_axis_off() -ax3.set_title('`match_template`\nresult') -# highlight matched region -ax3.autoscale(False) -ax3.plot(x, y, 'o', markeredgecolor='r', markerfacecolor='none', markersize=10) - -plt.show() - From 8a31b76ba8e5443d4adc2c7d53f053d38789f6b5 Mon Sep 17 00:00:00 2001 From: Tony S Yu Date: Tue, 8 May 2012 21:35:03 -0400 Subject: [PATCH 48/48] Add contribution note. --- CONTRIBUTORS.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 85043f2c..d4a3a816 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -28,7 +28,7 @@ - Tony Yu Reading of paletted images; build, bug and doc fixes. Code to generate skimage logo. - Otsu thresholding, histogram equalisation, and more. + Otsu thresholding, histogram equalisation, template matching, and more. - Zachary Pincus Tracing of low cost paths, FreeImage I/O plugin, iso-contours,