Python code to calculate color difference (read http://en.wikipedia.org/wiki/Color_difference for more information and the original formulae/algorithm) and convert RGB color space to Lab

coldiff.py

#!/usr/bin/env python

from math import sqrt


def rgb(x):
    """Convert #[AA]RRGGBB color in integer or string to (r,g,b) tuple
    
    Alpha (AA) component is simply ignored.
    
    rgb(0xff0000ff)
    >>> (0, 0, 255)
    rgb('#ff0000')
    >>> (255, 0, 0)
    """
    
    if isinstance(x, str) and x[0] == '#':
        x = int(x[1:], 16)
    return ((x >> 16) & 0xff, (x >> 8) & 0xff, (x) & 0xff)


def cie94(L1_a1_b1, L2_a2_b2):
    """Calculate color difference by using CIE94 formulae
    
    See http://en.wikipedia.org/wiki/Color_difference or
    http://www.brucelindbloom.com/index.html?Eqn_DeltaE_CIE94.html.
    
    cie94(rgb2lab((255, 255, 255)), rgb2lab((0, 0, 0)))
    >>> 58.0
    cie94(rgb2lab(rgb(0xff0000)), rgb2lab(rgb('#ff0000')))
    >>> 0.0
    """
    
    L1, a1, b1 = L1_a1_b1
    L2, a2, b2 = L2_a2_b2

    C1 = sqrt(_square(a1) + _square(b1))
    C2 = sqrt(_square(a2) + _square(b2))
    delta_L = L1 - L2
    delta_C = C1 - C2
    delta_a = a1 - a2
    delta_b = b1 - b2
    delta_H_square = _square(delta_a) + _square(delta_b) - _square(delta_C)
    return (sqrt(_square(delta_L)
            + _square(delta_C) / _square(1.0 + 0.045 * C1)
            + delta_H_square / _square(1.0 + 0.015 * C1)))


def rgb2lab(R_G_B):
    """Convert RGB colorspace to Lab
    
    Adapted from http://www.easyrgb.com/index.php?X=MATH.
    """
    
    R, G, B = R_G_B
    
    # Convert RGB to XYZ
    
    var_R = ( R / 255.0 )        # R from 0 to 255
    var_G = ( G / 255.0 )        # G from 0 to 255
    var_B = ( B / 255.0 )        # B from 0 to 255

    if ( var_R > 0.04045 ): var_R = ( ( var_R + 0.055 ) / 1.055 ) ** 2.4
    else:                   var_R = var_R / 12.92
    if ( var_G > 0.04045 ): var_G = ( ( var_G + 0.055 ) / 1.055 ) ** 2.4
    else:                   var_G = var_G / 12.92
    if ( var_B > 0.04045 ): var_B = ( ( var_B + 0.055 ) / 1.055 ) ** 2.4
    else:                   var_B = var_B / 12.92

    var_R = var_R * 100.0
    var_G = var_G * 100.0
    var_B = var_B * 100.0

    # Observer. = 2°, Illuminant = D65
    X = var_R * 0.4124 + var_G * 0.3576 + var_B * 0.1805
    Y = var_R * 0.2126 + var_G * 0.7152 + var_B * 0.0722
    Z = var_R * 0.0193 + var_G * 0.1192 + var_B * 0.9505
    
    # Convert XYZ to L*a*b*
    
    var_X = X / 95.047         # ref_X =  95.047   Observer= 2°, Illuminant= D65
    var_Y = Y / 100.000        # ref_Y = 100.000
    var_Z = Z / 108.883        # ref_Z = 108.883

    if ( var_X > 0.008856 ): var_X = var_X ** ( 1.0/3.0 )
    else:                    var_X = ( 7.787 * var_X ) + ( 16.0 / 116.0 )
    if ( var_Y > 0.008856 ): var_Y = var_Y ** ( 1.0/3.0 )
    else:                    var_Y = ( 7.787 * var_Y ) + ( 16.0 / 116.0 )
    if ( var_Z > 0.008856 ): var_Z = var_Z ** ( 1.0/3.0 )
    else:                    var_Z = ( 7.787 * var_Z ) + ( 16.0 / 116.0 )

    CIE_L = ( 116.0 * var_Y ) - 16.0
    CIE_a = 500.0 * ( var_X - var_Y )
    CIE_b = 200.0 * ( var_Y - var_Z )
    return (CIE_L, CIE_a, CIE_b)


def _old_rgb2lab(R_G_B):
    """Old implementation of rgb2lab, the result is strange :D
    
    Adapted from http://www.f4.fhtw-berlin.de/~barthel/ImageJ/ColorInspector//HTMLHelp/farbraumJava.htm.
    """
    
    R, G, B = R_G_B

    # http://www.brucelindbloom.com

    # float: r, g, b, X, Y, Z, fx, fy, fz, xr, yr, zr
    # float: Ls, as_, bs
    eps = 216.0/24389.0
    k = 24389.0/27.0

    Xr = 0.964221  # reference white D50
    Yr = 1.0
    Zr = 0.825211

    # RGB to XYZ
    r = R/255.0 #R 0..1
    g = G/255.0 #G 0..1
    b = B/255.0 #B 0..1
    
    if not (0<=r<=1 and 0<=g<=1 and 0<=b<=1):
        raise ValueError('RGB out of 0..255 range')

    # assuming sRGB (D65)
    if r <= 0.04045:
        r = r/12
    else:
        r = ((r+0.055)/1.055) ** 2.4

    if g <= 0.04045:
        g = g/12
    else:
        g = ((g+0.055)/1.055) ** 2.4

    if b <= 0.04045:
        b = b/12
    else:
        b = ((b+0.055)/1.055) ** 2.4

    X =  0.436052025*r     + 0.385081593*g + 0.143087414 *b
    Y =  0.222491598*r     + 0.71688606 *g + 0.060621486 *b
    Z =  0.013929122*r     + 0.097097002*g + 0.71418547  *b

    # XYZ to Lab
    xr = X/Xr
    yr = Y/Yr
    zr = Z/Zr

    if xr > eps:
        fx =  xr ** (1.0/3.0)
    else:
        fx = (k * xr + 16.0) / 116.0

    if yr > eps:
        fy =  yr ** (1.0/3.0)
    else:
        fy = (k * yr + 16.0) / 116.0

    if zr > eps:
        fz =  zr ** (1.0/3.0)
    else:
        fz = (k * zr + 16.0) / 116.0

    Ls = ( 116 * fy ) - 16
    as_ = 500*(fx-fy)
    bs = 200*(fy-fz)

    return (int(2.55*Ls + 0.5), # L
            int(as_ + 0.5),     # a
            int(bs + 0.5))      # b


def _square(x):
    return x * x