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':