diff --git a/skimage/data/lena_GRAY_U8_hog.npy b/skimage/data/lena_GRAY_U8_hog.npy
new file mode 100644
index 00000000..71ebf61f
Binary files /dev/null and b/skimage/data/lena_GRAY_U8_hog.npy differ
diff --git a/skimage/feature/_hog.py b/skimage/feature/_hog.py
index c0f4cba8..760102e7 100644
--- a/skimage/feature/_hog.py
+++ b/skimage/feature/_hog.py
@@ -5,7 +5,8 @@ from . import _hoghistogram
 
 
 def hog(image, orientations=9, pixels_per_cell=(8, 8),
-        cells_per_block=(3, 3), visualise=False, normalise=False):
+        cells_per_block=(3, 3), visualise=False, normalise=False,
+        feature_vector=True):
     """Extract Histogram of Oriented Gradients (HOG) for a given image.
 
     Compute a Histogram of Oriented Gradients (HOG) by
@@ -31,6 +32,9 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8),
     normalise : bool, optional
         Apply power law compression to normalise the image before
         processing.
+    feature_vector : bool, optional
+        Return the data as a feature vector by calling .ravel() on the result
+        just before returning.
 
     Returns
     -------
@@ -169,8 +173,11 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8),
     overlapping grid of blocks covering the detection window into a combined
     feature vector for use in the window classifier.
     """
+    
+    if feature_vector:
+        normalised_blocks = normalised_blocks.ravel()
 
     if visualise:
-        return normalised_blocks.ravel(), hog_image
+        return normalised_blocks, hog_image
     else:
-        return normalised_blocks.ravel()
+        return normalised_blocks
diff --git a/skimage/feature/tests/test_hog.py b/skimage/feature/tests/test_hog.py
index 6112d726..3afc0d89 100644
--- a/skimage/feature/tests/test_hog.py
+++ b/skimage/feature/tests/test_hog.py
@@ -1,5 +1,7 @@
+import os
 import numpy as np
 from scipy import ndimage as ndi
+import skimage as si
 from skimage import data
 from skimage import feature
 from skimage import img_as_float
@@ -9,7 +11,7 @@ from numpy.testing import (assert_raises,
                            )
 
 
-def test_histogram_of_oriented_gradients():
+def test_histogram_of_oriented_gradients_output_size():
     img = img_as_float(data.astronaut()[:256, :].mean(axis=2))
 
     fd = feature.hog(img, orientations=9, pixels_per_cell=(8, 8),
@@ -18,6 +20,17 @@ def test_histogram_of_oriented_gradients():
     assert len(fd) == 9 * (256 // 8) * (512 // 8)
 
 
+def test_histogram_of_oriented_gradients_output_correctness():
+    img = np.load(os.path.join(si.data_dir, 'lena_GRAY_U8.npy'))
+    correct_output = np.load(os.path.join(si.data_dir, 'lena_GRAY_U8_hog.npy'))
+    
+    output = feature.hog(img, orientations=9, pixels_per_cell=(8, 8), 
+                         cells_per_block=(3, 3), feature_vector=True,
+                         normalise=False, visualise=False)
+    
+    assert_almost_equal(output, correct_output)
+
+
 def test_hog_image_size_cell_size_mismatch():
     image = data.camera()[:150, :200]
     fd = feature.hog(image, orientations=9, pixels_per_cell=(8, 8),