Improved support for color images using lab transform.

skimage/exposure/exposure.py

@@ -2,7 +2,7 @@ import numpy as np
 
 import skimage
 from skimage.util.dtype import dtype_range
-from skimage.color import rgb2gray
+import skimage.color as color
 from skimage.util.dtype import convert
 
 from _adapthist import _adapthist
@@ -210,7 +210,7 @@ def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
 
     Returns
     -------
-    out - np.ndarray :
+    out : np.ndarray
         equalized image - may be a different shape than the original
 
     Notes
@@ -218,8 +218,17 @@ def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
     * The underlying algorithm relies on an image whose rows and columns are even multiples of
     the number of tiles, so the extra rows and columns are left at their original values, thus
     preserving the input image shape.
-    * For RGB or RGBA images, the algorithm is run on each channel.
+    * For grayscale images, CLAHE is performed on one channel, and a grayscale is returned
+    * For color images, the following steps are performed:
+       - The image is converted to LAB color space
+       - The CLAHE algorithm is run on the L channel
+       - The image is converted back to RGB space and returned
     * For RGBA images, the original alpha channel is removed.
+
+    References
+    ----------
+    .. [1] http://tog.acm.org/resources/GraphicsGems/
+    .. [2] https://en.wikipedia.org/wiki/CLAHE#CLAHE
     '''
     in_type = image.dtype.type
     if out_range == 'full':
@@ -227,23 +236,50 @@ def adapthist(image, nx=8, ny=8, clip_limit=0.01, nbins=256, out_range='full'):
     else:
         out_range = (image.min(), image.max())
     # must be converted to 12 bit for CLAHE
-    image = skimage.img_as_uint(image)
+    int_image = skimage.img_as_uint(image)
     MAX_VAL = 2 ** 12 - 1
-    image = rescale_intensity(image, out_range=(0, MAX_VAL))
+    int_image = rescale_intensity(int_image, out_range=(0, MAX_VAL))
     # handle color images - CLAHE accepts scalar images only
-    args = [image.copy(), 0, MAX_VAL, nx, ny, nbins, clip_limit]
+    args = [int_image.copy(), 0, MAX_VAL, nx, ny, nbins, clip_limit]
     if image.ndim == 3:
-        image = image[:, :, :3]
-        for channel in range(3):
-            args[0] = image[:, :, channel]
+        # check for grayscale
+        if (np.allclose(image[:, :, 0], image[:, :, 1]) and
+            np.allclose(image[:, :, 2], image[:, :, 3])):
+            args[0] = image[:, :, 0]
             out = _adapthist(*args)
-            image[:out.shape[0], :out.shape[1], channel] = out
+            image = int_image[:, :, :3]
+            for channel in range(3):
+                image[:out.shape[0], :out.shape[1], channel] = out
+        # for color images, convert to LAB space for processing
+        else:
+            lab_img = color.rgb2lab(skimage.img_as_float(image))
+            L_chan = lab_img[:, :, 0]
+            L_chan /= np.max(np.abs(L_chan))
+            L_chan = skimage.img_as_uint(L_chan)
+            args[0] = rescale_intensity(L_chan, out_range=(0, MAX_VAL))
+            new_L = _adapthist(*args).astype(float)
+            new_L = rescale_intensity(new_L, out_range=(0, 100))
+            lab_img[:new_L.shape[0], :new_L.shape[1], 0] = new_L
+            image = color.lab2rgb(lab_img)
+            image = rescale_intensity(image, out_range=(0, 1))
     else:
         out = _adapthist(*args)
+        image = int_image
         image[:out.shape[0], :out.shape[1]] = out
     # restore to desired output type and output limits
     image = rescale_intensity(image)
-    if in_type != np.uint16:
-        image = convert(image, in_type)
+    image = convert(image, in_type)
     image = rescale_intensity(image, out_range=out_range)
     return image
+
+if __name__ == '__main__':
+    from skimage import data
+    import matplotlib.pyplot as plt
+    img = skimage.img_as_uint(data.lena())
+    adapted = adapthist(img, nx=10, ny=9, clip_limit=0.01,
+                            nbins=128, out_range='original')
+    plt.imshow(img)
+    plt.figure(); plt.imshow(skimage.img_as_ubyte(adapted))
+    plt.figure(); plt.imshow(color.lab2rgb(color.rgb2lab(img)))
+    plt.show()
+    print 'Done'

skimage/exposure/tests/test_exposure.py

@@ -76,7 +76,7 @@ def test_rescale_out_range():
 # Test rescale intensity
 # ======================
 
-def test_adapthist_ubyte():
+def test_adapthist_scalar():
     '''Test a scalar uint8 image
     '''
     img = skimage.img_as_ubyte(data.moon())
@@ -84,23 +84,43 @@ def test_adapthist_ubyte():
     assert adapted.min() == 0
     assert adapted.max() == 255
     assert img.shape == adapted.shape
-    assert peak_snr(img, adapted) > 22
-    assert norm_brightness_err(img, adapted) < 0.05
+    full_scale = skimage.exposure.rescale_intensity(img)
+    assert_almost_equal = np.testing.assert_almost_equal
+    assert_almost_equal(peak_snr(full_scale, adapted), 22.19920073)
+    assert_almost_equal(norm_brightness_err(full_scale, adapted),
+                        0.04161278)
     return img, adapted
 
 
-def test_adapthist_float():
-    '''Test an RGB float image
+def test_adapthist_grayscale():
+    '''Test a grayscale float image
     '''
     img = skimage.img_as_float(data.lena())
+    img = rgb2gray(img)
     adapted = exposure.adapthist(img, nx=10, ny=9, clip_limit=0.01,
                         nbins=128, out_range='original')
     assert_almost_equal = np.testing.assert_almost_equal
     assert_almost_equal(adapted.min(), img.min())
-    assert_almost_equal(adapted.min(), img.min())
+    assert_almost_equal(adapted.max(), img.max())
     assert img.shape == adapted.shape
-    assert peak_snr(img, adapted) > 136
-    assert norm_brightness_err(img, adapted) < 0.02
+    assert_almost_equal(peak_snr(img, adapted), 131.4962063)
+    assert_almost_equal(norm_brightness_err(img, adapted), 0.0208805)
+    return data, adapted
+
+
+def test_adapthist_color():
+    '''Test a color uint16 image
+    '''
+    img = skimage.img_as_uint(data.lena())
+    adapted = exposure.adapthist(img, clip_limit=0.01)
+    assert_almost_equal = np.testing.assert_almost_equal
+    assert adapted.min() == 0
+    assert adapted.max() == 65535
+    assert img.shape == adapted.shape
+    full_scale = skimage.exposure.rescale_intensity(img)
+    assert_almost_equal(peak_snr(full_scale, adapted), 64.29546231)
+    assert_almost_equal(norm_brightness_err(full_scale, adapted),
+                        0.181473754)
     return data, adapted
 
 
@@ -135,7 +155,7 @@ def norm_brightness_err(img1, img2):
     Returns
     -------
     norm_brightness_error : float
-        Normalize absolute mean brightness error
+        Normalized absolute mean brightness error
     '''
     if img1.ndim == 3:
         img1, img2 = rgb2gray(img1), rgb2gray(img2)