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