diff --git a/skimage/measure/_marching_cubes_lewiner.py b/skimage/measure/_marching_cubes_lewiner.py
index 537395f2..f105b4b3 100644
--- a/skimage/measure/_marching_cubes_lewiner.py
+++ b/skimage/measure/_marching_cubes_lewiner.py
@@ -14,7 +14,8 @@ from . import _marching_cubes_lewiner_cy
 
 
 def marching_cubes_lewiner(volume, level=None, spacing=(1., 1., 1.),
-                           step_size=1, allow_degenerate=True, use_classic=False):
+                           gradient_direction='descent', step_size=1,
+                           allow_degenerate=True, use_classic=False):
     """
     Lewiner marching cubes algorithm to find surfaces in 3d volumetric data
     
@@ -31,14 +32,21 @@ def marching_cubes_lewiner(volume, level=None, spacing=(1., 1., 1.),
     spacing : length-3 tuple of floats
         Voxel spacing in spatial dimensions corresponding to numpy array
         indexing dimensions (M, N, P) as in `volume`.
+    gradient_direction : string
+        Controls if the mesh was generated from an isosurface with gradient
+        descent toward objects of interest (the default), or the opposite,
+        considering the *left-hand* rule.
+        The two options are:
+        * descent : Object was greater than exterior
+        * ascent : Exterior was greater than object
     step_size : int
         Step size in voxels. Default 1. Larger steps yield faster but
         coarser results. The result will always be topologically correct
         though.
     allow_degenerate : bool
-        Whether to allow degenerate triangles in the end-result. Default True.
-        If False, degenerate triangles are removed, making the algorithm 
-        about twice as slow.
+        Whether to allow degenerate (i.e. zero-area) triangles in the
+        end-result. Default True. If False, degenerate triangles are
+        removed, at the cost of making the algorithm slower.
     use_classic : bool
         If given and True, the classic marching cubes by Lorensen (1987)
         is used. This option is included for reference purposes. Note
@@ -120,6 +128,12 @@ def marching_cubes_lewiner(volume, level=None, spacing=(1., 1., 1.),
     
     # Finishing touches to output
     faces.shape = -1, 3
+    if gradient_direction == 'descent':
+        # MC implementation is right-handed, but gradient_direction is left-handed
+        faces = np.fliplr(faces)
+    elif not gradient_direction == 'ascent':
+        raise ValueError("Incorrect input %s in `gradient_direction`, see "
+                         "docstring." % (gradient_direction))
     if spacing != (1, 1, 1):
         vertices = vertices * np.r_[spacing]
     
diff --git a/skimage/measure/mc_meta/visual_test.py b/skimage/measure/mc_meta/visual_test.py
index ae8c0bfb..7a8ea699 100644
--- a/skimage/measure/mc_meta/visual_test.py
+++ b/skimage/measure/mc_meta/visual_test.py
@@ -9,10 +9,12 @@ import numpy as np
 import visvis as vv
 
 from skimage.measure import marching_cubes, marching_cubes_lewiner
+from skimage.draw import ellipsoid
 
 
 # Create test volume
-SELECT = 1
+SELECT = 4
+gradient_dir = 'descent'  # ascent or descent
 
 if SELECT == 1:
     # Medical data
@@ -45,20 +47,26 @@ elif SELECT == 3:
     # Uncommenting the line below will yield different results for classic MC
     #vol = -vol
 
+elif SELECT == 4:
+    vol = ellipsoid(4, 3, 2, levelset=True)
+    isovalue = 0
+
 # Get surface meshes
 t0 = time.time()
-vertices1, faces1, *_ = marching_cubes_lewiner(vol, isovalue, use_classic=False)
+vertices1, faces1, *_ = marching_cubes_lewiner(vol, isovalue, gradient_direction=gradient_dir, use_classic=False)
 print('finding surface lewiner took %1.0f ms' % (1000*(time.time()-t0)) )
 
 t0 = time.time()
-vertices2, faces2, *_ = marching_cubes(vol, isovalue)
+vertices2, faces2, *_ = marching_cubes(vol, isovalue, gradient_direction=gradient_dir)
 print('finding surface classic took %1.0f ms' % (1000*(time.time()-t0)) )
 
 # Show
 vv.figure(1); vv.clf()
-a1 = vv.subplot(121); vv.mesh(np.fliplr(vertices1), faces1)
-a2 = vv.subplot(122); vv.mesh(np.fliplr(vertices2), faces2)
+a1 = vv.subplot(121); m1 = vv.mesh(np.fliplr(vertices1), faces1)
+a2 = vv.subplot(122); m2 = vv.mesh(np.fliplr(vertices2), faces2)
 a1.camera = a2.camera
 
+# visvis uses right-hand rule, gradient_direction param uses left-hand rule
+m1.cullFaces = m2.cullFaces = 'front'  # None, front or back
 
 vv.use().Run()
diff --git a/skimage/measure/tests/test_marching_cubes.py b/skimage/measure/tests/test_marching_cubes.py
index 491fedfc..4d62abe0 100644
--- a/skimage/measure/tests/test_marching_cubes.py
+++ b/skimage/measure/tests/test_marching_cubes.py
@@ -106,21 +106,21 @@ def test_both_algs_same_result_ellipse():
     vertices2, faces2, *_ = marching_cubes_lewiner(sphere_small, 0, allow_degenerate=False)
     vertices3, faces3, *_ = marching_cubes_lewiner(sphere_small, 0, allow_degenerate=False, use_classic=True)
     
-    # Order id different, best we can do is test equal shape and same vertices present
+    # Order is different, best we can do is test equal shape and same vertices present
     assert _same_mesh(vertices1, faces1, vertices2, faces2)
     assert _same_mesh(vertices1, faces1, vertices3, faces3)
 
 
 def _same_mesh(vertices1, faces1, vertices2, faces2):
-    faces1 = [sorted(f) for f in faces1]
-    faces2 = [sorted(f) for f in faces2]
+    rounder = lambda x: int(x*1000)/1000  # to take into account small variations
     triangles1 = vertices1[np.array(faces1)]
     triangles2 = vertices2[np.array(faces2)]
-    triang1 = set([tuple(t.flat) for t in triangles1])
-    triang2 = set([tuple(t.flat) for t in triangles2])
+    triang1 = [np.concatenate(sorted(t, key=lambda x:tuple(x))) for t in triangles1]
+    triang1 = set([tuple([rounder(i) for i in t]) for t in triang1])
+    triang2 = [np.concatenate(sorted(t, key=lambda x:tuple(x))) for t in triangles2]
+    triang2 = set([tuple([rounder(i) for i in t]) for t in triang2])
     return triang1 == triang2
 
-
 def test_both_algs_same_result_donut():
     # Performing this test on data that does not have ambiguities