RGB-to-XYZ matrix calculation based on http://koichitamura.blogspot.com.au/2009/01/getting-xyz-tofrom-rgb-conversion.html

rgb_to_xyz_matrix.py

"""
RGB to XYZ gamut-remapping matrix calculation.

Main function is rgb_to_xyz_matrix.
"""

class MatrixError(Exception):
    pass

class NonInvertableMatrix(MatrixError):
    pass


class Matrix(object):
    def __init__(self, *matrix):
        if len(matrix) != 3*3:
            raise MatrixError(
                "Incorrect number of values, expected 3*3, got %d" % len(matrix))
        assert len(matrix) == 3*3
        self.matrix = matrix

    def __mul__(self, other):
        """Matrix'ify a vector3
        """
        r, g, b = other

        m = self.matrix
        nr = m[0] * r + m[1] * g + m[2] * b
        ng = m[3] * r + m[4] * g + m[5] * b
        nb = m[6] * r + m[7] * g + m[8] * b

        return (nr, ng, nb)

    def det(self):
        """Determinant of matrix
        http://en.wikipedia.org/wiki/Determinant#3-by-3_matrices
        """
        (a, b, c,
         d, e, f,
         g, h, i) = self.matrix

        thedet = a*e*i + b*f*g + c*d*h - a*f*h - b*d*i - c*e*g
        return thedet

    def invert(self):
        """Return an new, inverted matrix
        http://en.wikipedia.org/wiki/Matrix_inversion#Inversion_of_3.C3.973_matrices
        """

        (a, b, c,
         d, e, f,
         g, h, k) = self.matrix

        A = e*k - f*h
        B = f*g - d*k
        C = d*h - e*g
        D = c*h - b*k
        E = a*k - c*g
        F = b*g - a*h
        G = b*f - c*e
        H = c*d - a*f
        K = a*e - b*d

        new_matrix = [A, D, G, B, E, H, C, F, K]

        det = self.det()
        if det == 0:
            raise NonInvertableMatrix("Determinant was zero")

        new_matrix = [(1.0/det) * x for x in new_matrix]

        # Return new instance of self
        return self.__class__(*new_matrix)

    def __repr__(self):
        mvals = []
        for x in self.matrix:
            if x < 0:
                mvals.append("%.05f" % x)
            else:
                # Value has no -sign, so pad one space
                mvals.append(" %.05f" % x)

        outstr = "Matrix([\n"
        for x in range(3):
            outstr += "    " + ", ".join(mvals[x*3:x*3+3]) + ",\n"
        outstr += "])"
        return outstr


def assert_close(a, b, epsilon = 0.0001):
    """Test helper to verify float'ness, works on floats, or list/tuples of floats
    """
    if isinstance(a, (list, tuple)) and isinstance(b, (list, tuple)):
        for i in range(max(len(a), len(b))):
            assert_close(a[i], b[i])
    else:
        assert abs(a-b) < epsilon, "Difference between %r and %r > %r (difference is %r)" % (
            a, b, epsilon, abs(a-b))


def test_matrix_invert_uninvertable():
    """Tests we cannot invert an uninvertable matrix
    """
    m1 = Matrix(
        0, 0, 0,
        0, 1, 0,
        0, 0, 1)

    assert m1.det() == 0

    try:
        m2 = m1.invert()
    except NonInvertableMatrix:
        pass
    else:
        raise AssertionError("Expected NonInvertableMatrix exception, matrix should be non invertable, got: %s" % (repr(m2)))


def test_matrix_invert_identity():
    """Tests inverting an identity matrix gets an identity matrix
    """
    m1 = Matrix(
        1, 0, 0,
        0, 1, 0,
        0, 0, 1)

    assert m1.det() == 1

    m2 = m1.invert()

    print m1
    print m2

    assert m1.matrix == m2.matrix


def test_matrix_invert_almost_identity():
    """Tests inverting a simple matrix
    """
    m1 = Matrix(
        2, 0, 0,
        0, 1, 0,
        0, 0, 1)
    m2 = m1.invert()

    print m1
    print m2

    # First value should be 1/2
    assert_close(0.5, m2.matrix[0])

    # Rest of matrix should be same
    for a, b in zip(m1.matrix[1:], m2.matrix[1:]):
        assert_close(a, b)


def test_matrix_invert_more_complex():
    """Tests inverting a less trivial matrix
    """
    m1 = Matrix(
        2, 1, 0,
        0.1, 2.1, 1.1,
        0.1, 0.2, 2.2)

    print m1.det()
    assert_close(m1.det(), 8.69)

    m2 = m1.invert()
    assert_close(m2.det(), 0.115075)
    correct_matrix = (
         0.506329,  -0.253165,  0.126582,
        -0.0126582,  0.506329, -0.253165,
        -0.0218642, -0.0345224, 0.471807)

    for i, (a_b) in enumerate(zip(m2.matrix, correct_matrix)):
        a, b = a_b
        assert_close(a, b)


def rgb_to_xyz_matrix(red_xy, green_xy, blue_xy, refwhite_xyz, debug=False):
    """Based of this page:
    http://koichitamura.blogspot.com.au/2009/01/getting-xyz-tofrom-rgb-conversion.html

    Also somewhat useful, but either I misunderstood something or it contains an error somewhere:
    http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html
    """
    # Unpack primaries xy
    red_x, red_y = red_xy
    green_x, green_y = green_xy
    blue_x, blue_y = blue_xy

    # For each primary, get the "missing" z (where x+y+z = 1, so z = 1-x-y)
    red_z = 1 - red_xy[0] - red_xy[1]
    green_z = 1 - green_xy[0] - green_xy[1]rgb_to_xyz_matrix
    blue_z = 1 - blue_xy[0] - blue_xy[1]

    if debug:
        print("Red primary xyz: (%.04f, %0.4f, %.04f)" % (red_x, red_y, red_z))
        print("Green primary xyz: (%.04f, %0.4f, %.04f)" % (green_x, green_y, green_z))
        print("Blue primary xyz: (%.04f, %0.4f, %.04f)" % (blue_x, blue_y, blue_z))

    # Unpack refwhite XYZ
    refwhite_x, refwhite_y, refwhite_z = refwhite_xyz

    # Scale whitepoint xyz into XYZ, where Y equals 1.0
    refwhite_X = refwhite_x / refwhite_y
    refwhite_Y = refwhite_y / refwhite_y # redundant/same as Y=1
    refwhite_Z = refwhite_z / refwhite_y

    if debug:
        print("Reference whitepoint XYZ: (%.04f, %0.4f, %.04f)" % (refwhite_X, refwhite_Y, refwhite_Z))

    del refwhite_x, refwhite_y, refwhite_z # prevent stupid typos

    # Make a 3x3 matrix where the first row contains the red xyz
    # values etc
    K = Matrix(
        red_x, green_x, blue_x,
        red_y, green_y, blue_y,
        red_z, green_z, blue_z)

    if debug:
        print("K matrix: %r" % K)

    # Get scaling factors for the rgb primaries.  Done by inverting
    # the "K" matrix, and using it to matrixify the whitepoint XYZ
    scale_R, scale_G, scale_B = K.invert() * (refwhite_X, refwhite_Y, refwhite_Z)

    if debug:
        print("Scaling factors: %.04f, %.04f, %.04f" % (scale_R, scale_G, scale_B))

    # Make matrix with scaled primaries
    M = Matrix(
        red_x * scale_R, green_x * scale_G, blue_x * scale_B,
        red_y * scale_R, green_y * scale_G, blue_y * scale_B,
        red_z * scale_R, green_z * scale_G, blue_z * scale_B)

    # Done.
    return M


def test_srgb_to_xyz_matrix():
    """Values from:
    http://www.w3.org/Graphics/Color/sRGB

    sRGB primaries, D65 whitepoint

    Not quite the same as the ones in
    http://en.wikipedia.org/wiki/SRGB#The_sRGB_gamut
    ..which are allegedly in the sRGB spec
    """

    # Get the matrix to go from linearised sRGB values to XYZ
    m1 = rgb_to_xyz_matrix(
        red_xy = (0.64, 0.33),
        green_xy = (0.30, 0.60),
        blue_xy = (0.15, 0.06),
        refwhite_xyz = (0.3127, 0.3290, 0.3583))

    expected_srgb_to_xyz = (
        0.4124, 0.3576, 0.1805,
        0.2126, 0.7152, 0.0722,
        0.0193, 0.1192, 0.9505)

    for a, b in zip(m1.matrix, expected_srgb_to_xyz):
        assert_close(a, b)

    # Check the inverse matrix goes from XYZ to linearised sRGB
    m2 = m1.invert()
    expected_xyz_to_srgb = (
        3.241, -1.5374, -0.4986,
        -0.9692, 1.8760, 0.0416,
        0.0557, -0.2040, 1.0570)

    for a, b in zip(m2.matrix, expected_xyz_to_srgb):
        assert_close(a, b)


def test_matrix_logicalness():
    """Given an sRGB to XYZ matrix, matrixing rgb(1, 1, 1) will give
    the whitepoint, matrixing rgb(1, 0, 0) the red primary etc
    """

    refwhite = (0.3127, 0.3290, 0.3583)
    red_xy = (0.64, 0.33)
    green_xy = (0.30, 0.60)
    blue_xy = (0.15, 0.06)

    m1 = rgb_to_xyz_matrix(
        red_xy = red_xy,
        green_xy = green_xy,
        blue_xy = blue_xy,
        refwhite_xyz = refwhite)

    def XYZ_to_xyz(XYZ):
        # Scale XYZ so x+y+z = 1
        return [XYZ[i] / sum(XYZ) for i in range(len(XYZ))]

    outwhite = m1 * (1.0, 1.0, 1.0)
    outwhite = XYZ_to_xyz(outwhite)
    assert_close(1.0, sum(outwhite)) # Check the xyz values total 1
    assert_close(outwhite, refwhite) # and it matches the reference whitepoint

    red = m1 * (1.0, 0.0, 0.0)
    red = XYZ_to_xyz(red)
    assert_close(1.0, sum(red))
    assert_close(red[:2], red_xy) # Ignore z (it is verified by checking that x+y+z=1 on previous line)

    green = m1 * (0.0, 1.0, 0.0)
    green = XYZ_to_xyz(green)
    assert_close(1.0, sum(green))
    assert_close(green[:2], green_xy)

    blue = m1 * (0.0, 0.0, 1.0)
    blue = XYZ_to_xyz(blue)
    assert_close(1.0, sum(blue))
    assert_close(blue[:2], blue_xy)


if __name__ == "__main__":
    def diy_nosetest(testfunc):
        import sys
        import StringIO
        orig_stdout, sys.stdout = sys.stdout, StringIO.StringIO()

        try:
            testfunc()
        except Exception:
            sys.stdout, captured = orig_stdout, sys.stdout
            import traceback
            print traceback.format_exc()
            print "Captured stdout:"
            print captured.getvalue()
        finally:
            sys.stdout, _captured = orig_stdout, sys.stdout

    for x, func in locals().items():
        if x.startswith("test_"):
            diy_nosetest(func)


    print "\n\n# DCI P3"
    # Numbers from DCI specification 1.2

    dcip3_x = 0.3140
    dcip3_y = 0.3510

    # X = Y / y * Y
    # Y = 48 cd/m**2 (14 fL)
    # Z = Y / y * (1 - x - y)

    dcip3_Y = 48 # cd/m**2

    dcip3_X = (dcip3_Y / dcip3_y) * dcip3_x
    dcip3_Z = dcip3_Y / dcip3_y * (1 - dcip3_x - dcip3_y)

    dciP3_M = rgb_to_xyz_matrix(
        red_xy = (0.680, 0.320),
        green_xy = (0.265, 0.690),
        blue_xy = (0.150, 0.060),
        refwhite_xyz = (dcip3_X, dcip3_Y, dcip3_Z))

    print "DCI P3 to XYZ matrix:", dciP3_M
    print "XYZ to DCI P3 matrix:", dciP3_M.invert()


    # Some known primaries/whitepoints:
    # http://www.brucelindbloom.com/index.html?WorkingSpaceInfo.html
    # ..and the matricies:
    # http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html#WSMatrices

    print "\n\n# sRGB"

    srgb_M = rgb_to_xyz_matrix(
        red_xy = (0.64, 0.33),
        green_xy = (0.30, 0.60),
        blue_xy = (0.15, 0.06),
        refwhite_xyz = (0.3127, 0.3290, 0.3583))

    print "sRGB to XYZ matrix:", srgb_M
    print "XYZ to sRGB matrix:", srgb_M.invert()

    #white_XYZ = srgb_M * (1.0, 1.0, 1.0)

    #white_xyz = [white_XYZ[i] / sum(white_XYZ) for i in range(3)]
    #print white_xyz


    def XYZ_to_xy(X, Y, Z, with_z = False):
        x = X / sum((X, Y, Z))
        y = Y / sum((X, Y, Z))
        z = Z / sum((X, Y, Z))
        if with_z:
            return (x, y, z)
        else:
            return (x, y)

    # ACES to XYZ
    # From "ACES_1.0.1.pdf" https://github.com/ampas/aces-dev/tree/master/documents
    print "\n\n# ACES"
    aces_M = rgb_to_xyz_matrix(
        red_xy = (0.73470, 0.26530),
        green_xy = (0.0, 1.0),
        blue_xy = (0.00010, -0.07700),
        refwhite_xyz = (XYZ_to_xy(0.95265, 1.00000, 1.00883, with_z=True)))
    print "ACES to XYZ:", aces_M
    print "XYZ to ACES:", aces_M

    print "\n\n# Rec709"
    # http://en.wikipedia.org/wiki/Rec._709
    # with whitepoint XYZ from D65 page,
    # http://en.wikipedia.org/wiki/Illuminant_D65
    rec709_M = rgb_to_xyz_matrix(
        red_xy = (0.64, 0.33),
        green_xy = (0.30, 0.60),
        blue_xy = (0.15, 0.06),
        refwhite_xyz = (XYZ_to_xy(X=95.047, Y=100.00, Z=108.883, with_z=True))) # Illuminant D65
    print "Rec709 to XYZ", rec709_M
    print "XYZ to Rec709", rec709_M.invert()