diff --git a/doc/examples/plot_circular_hough_transform.py b/doc/examples/plot_circular_hough_transform.py
index fc891d04..d2f8f2ae 100755
--- a/doc/examples/plot_circular_hough_transform.py
+++ b/doc/examples/plot_circular_hough_transform.py
@@ -66,7 +66,7 @@ for idx in np.argsort(accums)[::-1][:5]:
     center_x, center_y = centers[idx]
     radius = radii[idx]
     cx, cy = circle_perimeter(center_y, center_x, radius)
-    image[cy, cx] = (220, 250, 150)
+    image[cy, cx] = (220, 20, 20)
 
 ax.imshow(image, cmap=plt.cm.gray)
 plt.show()
diff --git a/skimage/transform/_hough_transform.pyx b/skimage/transform/_hough_transform.pyx
index 29a38eb7..2b1acc7c 100644
--- a/skimage/transform/_hough_transform.pyx
+++ b/skimage/transform/_hough_transform.pyx
@@ -20,10 +20,9 @@ cdef inline Py_ssize_t round(double r):
     return <Py_ssize_t>((r + 0.5) if (r > 0.0) else (r - 0.5))
 
 
-@cython.boundscheck(False)
-def _hough_circle(cnp.ndarray img, \
-                  cnp.ndarray[ndim=1, dtype=cnp.intp_t] radius, \
-                  normalize=True):
+def _hough_circle(cnp.ndarray img,
+                  cnp.ndarray[ndim=1, dtype=cnp.intp_t] radius,
+                  char normalize=True):
     """Perform a circular Hough transform.
 
     Parameters
@@ -49,32 +48,42 @@ def _hough_circle(cnp.ndarray img, \
     cdef cnp.ndarray[ndim=1, dtype=cnp.intp_t] x, y
     x, y = np.nonzero(img)
 
+    cdef Py_ssize_t num_pixels = x.size
+
     # Offset the image
-    cdef int max_radius = radius.max()
+    cdef Py_ssize_t max_radius = radius.max()
     x = x + max_radius
     y = y + max_radius
 
-    cdef int px, py
-    cdef cnp.ndarray[ndim=1, dtype=cnp.intp_t] tx, ty, circle_x, circle_y
-    cdef cnp.ndarray acc = np.zeros((radius.size,
-                                     img.shape[0] + 2 * max_radius,
-                                     img.shape[1] + 2 * max_radius))
+    cdef Py_ssize_t i, p, c, num_circle_pixels, tx, ty
+    cdef double incr
+    cdef cnp.ndarray[ndim=1, dtype=cnp.intp_t] circle_x, circle_y
+
+    cdef cnp.ndarray[ndim=3, dtype=cnp.double_t] acc = \
+         np.zeros((radius.size,
+                   img.shape[0] + 2 * max_radius,
+                   img.shape[1] + 2 * max_radius), dtype=np.double)
 
     for i, rad in enumerate(radius):
         # Store in memory the circle of given radius
         # centered at (0,0)
         circle_x, circle_y = circle_perimeter(0, 0, rad)
 
-        # For each non zero pixel
-        for (px, py) in zip(x, y):
-            # Plug the circle at (px, py),
-            # its coordinates are (tx, ty)
-            tx = circle_x + px
-            ty = circle_y + py
-            acc[i, tx, ty] += 1
+        num_circle_pixels = circle_x.size
 
         if normalize:
-            acc[i] = acc[i] / len(circle_x)
+            incr = 1.0 / num_circle_pixels
+        else:
+            incr = 1
+
+        # For each non zero pixel
+        for p in range(num_pixels):
+            # Plug the circle at (px, py),
+            # its coordinates are (tx, ty)
+            for c in range(num_circle_pixels):
+                tx = circle_x[c] + x[p]
+                ty = circle_y[c] + y[p]
+                acc[i, tx, ty] += incr
 
     return acc