Create LAB LUTs for skimage.color.colorconv based on OpenCV Method

lab_tables.py

import numpy as np
from scipy.interpolate import interp1d
from skimage import img_as_ubyte


def build_spline(x, y):
    """computes cubic spline coefficients for a function: (xi=i, yi=f[i]), i=0..n

    Adapted rom OpenCV project.
    """
    func = interp1d(x, y, kind='cubic')
    xnew = np.linspace(0, 1, x.size * 4)
    return func(xnew)


LAB_CBRT_TAB_SIZE = 1024
GAMMA_TAB_SIZE = 1024
LabCbrtTabScale = LAB_CBRT_TAB_SIZE / 1.5
GammaTabScale = GAMMA_TAB_SIZE
XYZ_SHIFT = 12
LAB_SHIFT = XYZ_SHIFT
GAMMA_SHIFT = 3
LAB_SHIFT2 = LAB_SHIFT + GAMMA_SHIFT
LAB_CBRT_TAB_SIZE_B = (256 * 3 / 2. * (1 << GAMMA_SHIFT))
sRGB2XYZ_D65 = [0.412453, 0.357580, 0.180423,
                0.212671, 0.715160, 0.072169,
                0.019334, 0.119193, 0.950227]
D65 = [0.950456, 1., 1.088754]


def create_lab_tables():
    """Compute lab LUTs

    Adapted from OpenCV Project.
    """
    f = np.zeros(LAB_CBRT_TAB_SIZE + 1)
    g = np.zeros(GAMMA_TAB_SIZE + 1)
    ig = np.zeros(GAMMA_TAB_SIZE + 1)

    x = np.linspace(0, 1, f.size) * LabCbrtTabScale / LAB_CBRT_TAB_SIZE
    mask = x < 0.008856
    f[mask] = x[mask] * 7.787 + 0.13793103448275862
    f[~mask] = np.power(x[~mask], 1 / 3.)

    lab_cbrt_tab = x  #build_spline(x, f)

    x = np.linspace(0, 1, g.size)
    mask = x <= 0.04045
    g[mask] = x[mask] / 12.92
    g[~mask] = np.power((x[~mask] + 0.055) / 1.055, 2.4)

    mask = x <= 0.0031308
    ig[mask] = x[mask] * 12.92
    ig[~mask] = 1.055 * np.power(x[~mask], 1 / 2.4) - 0.055

    sRGB_gamma_tab = x  #build_spline(x, g)
    sRGB_inv_gamma_tab = x  # build_spline(x, ig)

    x = np.linspace(0, 1, 256)
    mask = x <= 0.04045

    sRGB_gamma_tab_b = np.ones(256) * 255 * (1 << GAMMA_SHIFT)
    sRGB_gamma_tab_b[mask] *= x[mask] / 12.92
    sRGB_gamma_tab_b[~mask] *= np.power((x[~mask] + 0.055) / 1.055, 2.4)
    sRGB_gamma_tab_b[sRGB_gamma_tab_b > 2 ** 16 - 1] = 2 ** 16 - 1
    sRGB_gamma_tab_b = sRGB_gamma_tab_b.astype(np.int32)

    linear_gamma_tab_b = np.arange(256, dtype=np.int32) * (1 << GAMMA_SHIFT)

    x = np.arange(LAB_CBRT_TAB_SIZE_B) / (255. * (1 << GAMMA_SHIFT))
    mask = x < 0.008856

    lab_cbrt_tab_b = np.zeros(LAB_CBRT_TAB_SIZE_B, dtype=np.int32)
    lab_cbrt_tab_b[mask] = ((1 << LAB_SHIFT2) *
        (x[mask] * 7.787 + 0.13793103448275862))
    lab_cbrt_tab_b[~mask] = (1 << LAB_SHIFT2) * np.power(x[~mask], 1 / 3.)
    lab_cbrt_tab_b[lab_cbrt_tab_b > 2 ** 16 - 1] = 2 ** 16 - 1
    lab_cbrt_tab_b = lab_cbrt_tab_b.astype(np.int32)

    return (sRGB_gamma_tab, sRGB_inv_gamma_tab, sRGB_gamma_tab_b,
        linear_gamma_tab_b, lab_cbrt_tab, lab_cbrt_tab_b)

(sRGB_gamma_tab, sRGB_inv_gamma_tab, sRGB_gamma_tab_b,
        linear_gamma_tab_b, lab_cbrt_tab, lab_cbrt_tab_b) = create_lab_tables()


@profile
def rgb2lab(image, coeffs=None, whitept=None):

    image = img_as_ubyte(image)

    blueIdx = 2  # for RGB

    if coeffs is None:
        coeffs = sRGB2XYZ_D65
    if whitept is None:
        whitept = D65
    scale = [float(1 << LAB_SHIFT) / whitept[0],
                     float(1 << LAB_SHIFT),
                    float(1 << LAB_SHIFT) / whitept[2]]

    C = np.zeros(9, dtype=np.int32)
    for i in range(3):
        j = i * 3
        C[j + np.bitwise_xor(blueIdx, 2)] = np.round(coeffs[j] * scale[i])
        C[j + 1] = np.round(coeffs[j + 1] * scale[i])
        C[j + blueIdx] = np.round(coeffs[j + 2] * scale[i])
        assert (C[j] >= 0 and C[j + 1] >= 0 and C[j + 2] >= 0
            and C[j] + C[j + 1] + C[j + 2] < 2 * (1 << LAB_SHIFT))

    Lscale = (116*255+50)//100
    Lshift = -((16*255*(1 << LAB_SHIFT2) + 50) // 100)

    tab = sRGB_gamma_tab_b

    def descale(x, n):
        return (((x) + (1 << ((n)-1))) >> (n))

    data = np.take(tab, image)
    dot = np.einsum('ijk,...k', data, C.reshape(3, 3))

    XYZ = descale(dot, LAB_SHIFT)

    fX = lab_cbrt_tab_b[XYZ[0]]
    fY = lab_cbrt_tab_b[XYZ[1]]
    fZ = lab_cbrt_tab_b[XYZ[2]]

    L = descale((Lscale * fY + Lshift), LAB_SHIFT2)
    L[L > 255] = 255
    offset = 128 * (1 << LAB_SHIFT2)
    a = descale((500 * (fX - fY) + offset), LAB_SHIFT2)
    a[a > 255] = 255
    b = descale((200 * (fY - fZ) + offset), LAB_SHIFT2)
    b[b > 255] = 255
    image[:, :, 0] = L
    image[:, :, 1] = a
    image[:, :, 2] = b

    return image.astype(np.uint8)


import time
from skimage import data, color
import cv2

image = data.chelsea()

t = time.time()
lab = color.rgb2lab(image.copy())
print time.time() - t

t = time.time()
lab2 = rgb2lab(image.copy())
print time.time() - t

t = time.time()
lab3 = cv2.cvtColor(image.copy(), cv2.cv.CV_RGB2Lab)
print time.time() - t

rgb2 = cv2.cvtColor(lab2, cv2.cv.CV_Lab2RGB)

np.testing.assert_allclose(lab2, lab3, rtol=1)
assert np.sum(rgb2.astype(float) - image) / image.size < 0.015

import matplotlib.pyplot as plt
plt.subplot(211)
plt.imshow(image)
plt.subplot(212)
plt.imshow(rgb2)
plt.tight_layout()
plt.show()