Move template matching to feature subpackage

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

skimage/feature/_template.pyx

@@ -0,0 +1,121 @@
+"""template.py - Template matching
+"""
+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":
+    double sqrt(double x)
+    double fabs(double x)
+
+
+@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.
+
+    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
+        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[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),
+                                              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.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])
+    # 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
+

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
 

skimage/feature/template.py

@@ -0,0 +1,84 @@
+"""template.py - Template matching
+"""
+import numpy as np
+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"):
+    """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.
+
+    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)
+

skimage/feature/tests/test_template.py

@@ -0,0 +1,47 @@
+import numpy as np
+from skimage.feature import match_template
+from numpy.random import randn
+
+
+def test_template():
+    size = 100
+    image = np.zeros((400, 400))
+    target = np.tri(size) + np.tri(size)[::-1]
+    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__":
+    from numpy import testing
+    testing.run_module_suite()
+