diff --git a/skimage/measure/_moments.pyx b/skimage/measure/_moments.pyx
index e6ccdb75..5a9bc11a 100644
--- a/skimage/measure/_moments.pyx
+++ b/skimage/measure/_moments.pyx
@@ -5,7 +5,7 @@
 import numpy as np
 
 
-def moments(double[:, :] image, Py_ssize_t order=3):
+cpdef moments(double[:, :] image, Py_ssize_t order=3):
     """Calculate all raw image moments up to a certain order.
 
     The following properties can be calculated from raw image moments:
@@ -42,8 +42,8 @@ def moments(double[:, :] image, Py_ssize_t order=3):
     return moments_central(image, 0, 0, order)
 
 
-def moments_central(double[:, :] image, double cr, double cc,
-                    Py_ssize_t order=3):
+cpdef moments_central(double[:, :] image, double cr, double cc,
+                      Py_ssize_t order=3):
     """Calculate all central image moments up to a certain order.
 
     Note that central moments are translation invariant but not scale and
@@ -79,11 +79,19 @@ def moments_central(double[:, :] image, double cr, double cc,
     """
     cdef Py_ssize_t p, q, r, c
     cdef double[:, ::1] mu = np.zeros((order + 1, order + 1), dtype=np.double)
-    for p in range(order + 1):
-        for q in range(order + 1):
-            for r in range(image.shape[0]):
-                for c in range(image.shape[1]):
-                    mu[p, q] += image[r, c] * (r - cr) ** q * (c - cc) ** p
+    cdef double val, dr, dc, dcp, drq
+    for r in range(image.shape[0]):
+        dr = r - cr
+        for c in range(image.shape[1]):
+            dc = c - cc
+            val = image[r, c]
+            dcp = 1
+            for p in range(order + 1):
+                drq = 1
+                for q in range(order + 1):
+                    mu[p, q] += val * drq * dcp
+                    drq *= dr
+                dcp *= dc
     return np.asarray(mu)