diff --git a/skimage/feature/_hog.py b/skimage/feature/_hog.py
index 11f7099d..b9d62926 100644
--- a/skimage/feature/_hog.py
+++ b/skimage/feature/_hog.py
@@ -2,7 +2,7 @@ import numpy as np
 from scipy import sqrt, pi, arctan2, cos, sin
 from scipy.ndimage import uniform_filter
 from .._shared.utils import assert_nD
-
+import _hoghistogram
 
 def hog(image, orientations=9, pixels_per_cell=(8, 8),
         cells_per_block=(3, 3), visualise=False, normalise=False):
@@ -114,24 +114,10 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8),
     n_cellsx = int(np.floor(sx // cx))  # number of cells in x
     n_cellsy = int(np.floor(sy // cy))  # number of cells in y
 
-    # compute orientations integral images
     orientation_histogram = np.zeros((n_cellsy, n_cellsx, orientations))
-    subsample = np.index_exp[cy // 2:cy * n_cellsy:cy,
-                             cx // 2:cx * n_cellsx:cx]
-    for i in range(orientations):
-        # create new integral image for this orientation
-        # isolate orientations in this range
 
-        temp_ori = np.where(orientation < 180.0 / orientations * (i + 1),
-                            orientation, -1)
-        temp_ori = np.where(orientation >= 180.0 / orientations * i,
-                            temp_ori, -1)
-        # select magnitudes for those orientations
-        cond2 = temp_ori > -1
-        temp_mag = np.where(cond2, magnitude, 0)
-
-        temp_filt = uniform_filter(temp_mag, size=(cy, cx))
-        orientation_histogram[:, :, i] = temp_filt[subsample]
+    _hoghistogram.HogHistograms(gx, gy, cx, cy, sx, sy, n_cellsx, n_cellsy, visualise, orientations, 
+        orientation_histogram)
 
     # now for each cell, compute the histogram
     hog_image = None
diff --git a/skimage/feature/_hoghistogram.pyx b/skimage/feature/_hoghistogram.pyx
new file mode 100644
index 00000000..96112730
--- /dev/null
+++ b/skimage/feature/_hoghistogram.pyx
@@ -0,0 +1,70 @@
+# cython: profile=True
+# cython: cdivision=True
+# cython: boundscheck=False
+# cython: wraparound=False
+
+import cmath, math
+cimport numpy as np
+import numpy as np
+from scipy import pi, arctan2, cos, sin
+
+cdef float CellHog(np.ndarray[np.float64_t, ndim=2] magnitude, 
+	np.ndarray[np.float64_t, ndim=2] orientation,
+	float ori1, float ori2,
+	int cx, int cy, int xi, int yi, int sx, int sy):
+	cdef int cx1, cy1
+
+	cdef float total = 0.
+	for cy1 in range(-cy/2, cy/2):
+		for cx1 in range(-cx/2, cx/2):
+			if yi + cy1 < 0: continue
+			if yi + cy1 >= sy: continue
+			if xi + cx1 < 0: continue
+			if xi + cx1 >= sx: continue
+			if orientation[yi + cy1, xi + cx1] >= ori1: continue
+			if orientation[yi + cy1, xi + cx1] < ori2: continue
+
+			total += magnitude[yi + cy1, xi + cx1]
+
+	return total
+
+def HogHistograms(np.ndarray[np.float64_t, ndim=2] gx, \
+	np.ndarray[np.float64_t, ndim=2] gy, 
+	int cx, int cy, #Pixels per cell
+	int sx, int sy, #Image size
+	int n_cellsx, int n_cellsy, 
+	int visualise, int orientations, 
+	np.ndarray[np.float64_t, ndim=3] orientation_histogram):
+
+	cdef np.ndarray[np.float64_t, ndim=2] magnitude = np.sqrt(gx**2 + gy**2)
+	cdef np.ndarray[np.float64_t, ndim=2] orientation = arctan2(gy, gx) * (180 / pi) % 180
+	cdef int i, x, y, o, yi, xi, cy1, cy2, cx1, cx2
+	cdef float ori1, ori2
+
+	# compute orientations integral images
+
+	for i in range(orientations):
+		# isolate orientations in this range
+
+		ori1 = 180. / orientations * (i + 1)
+		ori2 = 180. / orientations * i
+
+		y = cy / 2
+		cy2 = cy * n_cellsy
+		x = cx / 2
+		cx2 = cx * n_cellsx
+		yi = 0
+		xi = 0
+
+		while y < cy2:
+			xi = 0
+			x = cx / 2
+
+			while x < cx2:
+				orientation_histogram[yi, xi, i] = CellHog(magnitude, orientation, ori1, ori2, cx, cy, x, y, sx, sy)
+				xi += 1
+				x += cx
+
+			yi += 1
+			y += cy
+
diff --git a/skimage/feature/setup.py b/skimage/feature/setup.py
index d7be3dcf..62bd950b 100644
--- a/skimage/feature/setup.py
+++ b/skimage/feature/setup.py
@@ -18,6 +18,7 @@ def configuration(parent_package='', top_path=None):
     cython(['brief_cy.pyx'], working_path=base_path)
     cython(['_texture.pyx'], working_path=base_path)
     cython(['_hessian_det_appx.pyx'], working_path=base_path)
+    cython(['_hoghistogram.pyx'], working_path=base_path)
 
     config.add_extension('corner_cy', sources=['corner_cy.c'],
                          include_dirs=[get_numpy_include_dirs()])
@@ -31,6 +32,8 @@ def configuration(parent_package='', top_path=None):
                          include_dirs=[get_numpy_include_dirs(), '../_shared'])
     config.add_extension('_hessian_det_appx', sources=['_hessian_det_appx.c'],
                          include_dirs=[get_numpy_include_dirs()])
+    config.add_extension('_hoghistogram', sources=['_hoghistogram.c'],
+                         include_dirs=[get_numpy_include_dirs(), '../_shared'])
 
     return config