diff --git a/skimage/color/colorconv.py b/skimage/color/colorconv.py
index 8e8c3b41..938bb849 100644
--- a/skimage/color/colorconv.py
+++ b/skimage/color/colorconv.py
@@ -53,6 +53,7 @@ References
 
 from __future__ import division
 
+from warnings import warn
 import numpy as np
 from scipy import linalg
 from ..util import dtype
@@ -552,6 +553,8 @@ def xyz2rgb(xyz):
     mask = arr > 0.0031308
     arr[mask] = 1.055 * np.power(arr[mask], 1 / 2.4) - 0.055
     arr[~mask] *= 12.92
+    arr[arr < 0] = 0
+    arr[arr > 1] = 1
     return arr
 
 
@@ -808,7 +811,7 @@ def xyz2lab(xyz, illuminant="D65", observer="2"):
     return np.concatenate([x[..., np.newaxis] for x in [L, a, b]], axis=-1)
 
 
-def lab2xyz(lab, illuminant="D65", observer="2"):
+def lab2xyz(lab, illuminant="D65", observer="2", error_on_invalid=False):
     """CIE-LAB to XYZcolor space conversion.
 
     Parameters
@@ -819,6 +822,9 @@ def lab2xyz(lab, illuminant="D65", observer="2"):
         The name of the illuminant (the function is NOT case sensitive).
     observer : {"2", "10"}, optional
         The aperture angle of the observer.
+     error_on_invalid: optional, bool
+        If true, raise ValueError when pixels are outside the valid gamut.
+        Otherwise, raise a warning.
 
     Returns
     -------
@@ -854,6 +860,15 @@ def lab2xyz(lab, illuminant="D65", observer="2"):
     x = (a / 500.) + y
     z = y - (b / 200.)
 
+    if np.any(z < 0):
+        invalid = np.nonzero(z < 0)
+        msg = 'Color data out of range: Z < 0 in %s pixels' % invalid[0].size
+        if error_on_invalid:
+            raise ValueError(msg)
+        else:
+            warn(msg)
+            z[invalid] = 0
+
     out = np.dstack([x, y, z])
 
     mask = out > 0.2068966
diff --git a/skimage/color/tests/test_colorconv.py b/skimage/color/tests/test_colorconv.py
index 094873ae..2b435f67 100644
--- a/skimage/color/tests/test_colorconv.py
+++ b/skimage/color/tests/test_colorconv.py
@@ -11,6 +11,7 @@ Authors
 :license: modified BSD
 """
 
+from __future__ import division
 import os.path
 
 import numpy as np
@@ -372,6 +373,29 @@ class TestColorconv(TestCase):
         img_rgb = img_as_float(self.img_rgb)
         assert_array_almost_equal(luv2rgb(rgb2luv(img_rgb)), img_rgb)
 
+    def test_lab_rgb_outlier(self):
+        lab_array = np.ones((3, 1, 3))
+        lab_array[0] = [50, -12, 85]
+        lab_array[1] = [50, 12, -85]
+        lab_array[2] = [90, -4, -47]
+        rgb_array = np.array([[[0.501, 0.481, 0]],
+                              [[0, 0.482, 1.]],
+                              [[0.578, 0.914, 1.]],
+                              ])
+        assert_almost_equal(lab2rgb(lab_array), rgb_array, decimal=3)
+
+    def test_lab_full_gamut(self):
+        a, b = np.meshgrid(np.arange(-100, 100), np.arange(-100, 100))
+        L = np.ones(a.shape)
+        lab = np.dstack((L, a, b))
+        assert_raises(ValueError, lambda: lab2xyz(lab))
+        multipliers = [0, 10, 50, 100]
+        nan_percents = [12, 9, 0, 0]
+        for (l, perc) in zip(multipliers, nan_percents):
+            lab[:, :, 0] = l
+            out = lab2xyz(lab)
+            assert_equal(int(np.isnan(out).sum() / out.size * 100), perc)
+
     def test_lab_lch_roundtrip(self):
         rgb = img_as_float(self.img_rgb)
         lab = rgb2lab(rgb)