Express colour conversions as chromaticities and adjust so that

everything is specified as something_to_xyz and then you can get
an inverse LUT if you want one.

16 files changed, 377 insertions, 122 deletions

diff --git a/examples/read_dcp.cc b/examples/read_dcp.cc
index d1416b28..565e6e4e 100644
--- a/examples/read_dcp.cc
+++ b/examples/read_dcp.cc
@@ -88,7 +88,7 @@ main ()
 
 	/* Convert to ARGB */
 	boost::scoped_array<uint8_t> rgba (new uint8_t[picture_image_xyz->size().width * picture_image_xyz->size().height * 4]);
-	dcp::xyz_to_rgba (picture_image_xyz, dcp::ColourConversion::xyz_to_srgb(), rgba.get ());
+	dcp::xyz_to_rgba (picture_image_xyz, dcp::ColourConversion::srgb_to_xyz(), rgba.get ());
 
 	Magick::Image image (picture_image_xyz->size().width, picture_image_xyz->size().height, "BGRA", Magick::CharPixel, rgba.get ());
 	image.write ("frame.png");
diff --git a/src/chromaticity.h b/src/chromaticity.h
new file mode 100644
index 00000000..89b71936
--- /dev/null
+++ b/src/chromaticity.h
@@ -0,0 +1,54 @@
+/*
+    Copyright (C) 2015 Carl Hetherington <cth@carlh.net>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program; if not, write to the Free Software
+    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+
+*/
+
+#ifndef DCP_CHROMATICITY_H
+#define DCP_CHROMATICITY_H
+
+#include <cmath>
+
+namespace dcp {
+
+class Chromaticity
+{
+public:
+	Chromaticity ()
+		: x (0)
+		, y (0)
+	{}
+
+	Chromaticity (double x_, double y_)
+		: x (x_)
+		, y (y_)
+	{}
+		
+	double x;
+	double y;
+
+	double z () const {
+		return 1 - x - y;
+	}
+
+	bool about_equal (Chromaticity const & other, float epsilon) const {
+		return std::fabs (x - other.x) < epsilon && std::fabs (y - other.y) < epsilon;
+	}
+};
+
+}
+
+#endif
diff --git a/src/colour_conversion.cc b/src/colour_conversion.cc
index 39799927..a5f77f0a 100644
--- a/src/colour_conversion.cc
+++ b/src/colour_conversion.cc
@@ -21,28 +21,45 @@
 #include "gamma_transfer_function.h"
 #include "modified_gamma_transfer_function.h"
 #include "colour_matrix.h"
+#include "dcp_assert.h"
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/lu.hpp>
+#include <boost/numeric/ublas/io.hpp>
 
 using boost::shared_ptr;
+using boost::optional;
 using namespace dcp;
 
 ColourConversion const &
 ColourConversion::srgb_to_xyz ()
 {
 	static ColourConversion* c = new ColourConversion (
-		shared_ptr<const TransferFunction> (new ModifiedGammaTransferFunction (false, 2.4, 0.04045, 0.055, 12.92)),
-		dcp::colour_matrix::rgb_to_xyz,
-		shared_ptr<const TransferFunction> (new GammaTransferFunction (true, 2.6))
+		shared_ptr<const TransferFunction> (new ModifiedGammaTransferFunction (2.4, 0.04045, 0.055, 12.92)),
+		YUV_TO_RGB_REC601,
+		Chromaticity (0.64, 0.33),
+		Chromaticity (0.3, 0.6),
+		Chromaticity (0.15, 0.06),
+		/* D65 */
+		Chromaticity (0.3127, 0.329),
+		optional<Chromaticity> (),
+		shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6))
 		);
 	return *c;
 }
 
 ColourConversion const &
-ColourConversion::xyz_to_srgb ()
+ColourConversion::rec601_to_xyz ()
 {
 	static ColourConversion* c = new ColourConversion (
-		shared_ptr<const TransferFunction> (new GammaTransferFunction (false, 2.6)),
-		dcp::colour_matrix::xyz_to_rgb,
-		shared_ptr<const TransferFunction> (new ModifiedGammaTransferFunction (true, 2.4, 0.04045, 0.055, 12.92))
+		shared_ptr<const TransferFunction> (new ModifiedGammaTransferFunction (1 / 0.45, 0.081, 0.099, 4.5)),
+		YUV_TO_RGB_REC601,
+		Chromaticity (0.64, 0.33),
+		Chromaticity (0.3, 0.6),
+		Chromaticity (0.15, 0.06),
+		/* D65 */
+		Chromaticity (0.3127, 0.329),
+		optional<Chromaticity> (),
+		shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6))
 		);
 	return *c;
 }
@@ -51,39 +68,169 @@ ColourConversion const &
 ColourConversion::rec709_to_xyz ()
 {
 	static ColourConversion* c = new ColourConversion (
-		shared_ptr<const TransferFunction> (new ModifiedGammaTransferFunction (false, 1 / 0.45, 0.081, 0.099, 4.5)),
-		dcp::colour_matrix::rgb_to_xyz,
-		shared_ptr<const TransferFunction> (new GammaTransferFunction (true, 2.6))
+		shared_ptr<const TransferFunction> (new ModifiedGammaTransferFunction (1 / 0.45, 0.081, 0.099, 4.5)),
+		YUV_TO_RGB_REC709,
+		Chromaticity (0.64, 0.33),
+		Chromaticity (0.3, 0.6),
+		Chromaticity (0.15, 0.06),
+		/* D65 */
+		Chromaticity (0.3127, 0.329),
+		optional<Chromaticity> (),
+		shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6))
 		);
 	return *c;
 }
 
 ColourConversion::ColourConversion (
 	shared_ptr<const TransferFunction> in,
-	double const matrix[3][3],
+	YUVToRGB yuv_to_rgb,
+	Chromaticity red,
+	Chromaticity green,
+	Chromaticity blue,
+	Chromaticity white,
+	optional<Chromaticity> adjusted_white,
 	shared_ptr<const TransferFunction> out
 	)
 	: _in (in)
-	, _matrix (3, 3)
+	, _yuv_to_rgb (yuv_to_rgb)
+	, _red (red)
+	, _green (green)
+	, _blue (blue)
+	, _white (white)
+	, _adjusted_white (adjusted_white)
 	, _out (out)
 {
-	for (int i = 0; i < 3; ++i) {
-		for (int j = 0; j < 3; ++j) {
-			_matrix (i, j) = matrix[i][j];
-		}
-	}
+	
 }
 
 bool
 ColourConversion::about_equal (ColourConversion const & other, float epsilon) const
 {
-	for (int i = 0; i < 3; ++i) {
-		for (int j = 0; j < 3; ++j) {
-			if (fabs (_matrix(i, j) - other._matrix(i, j)) > epsilon) {
-				return false;
-			}
-		}
+	if (!_in->about_equal (other._in, epsilon) ||
+	    _yuv_to_rgb != other._yuv_to_rgb ||
+	    !_red.about_equal (other._red, epsilon) ||
+	    !_green.about_equal (other._green, epsilon) ||
+	    !_blue.about_equal (other._blue, epsilon) ||
+	    !_white.about_equal (other._white, epsilon) ||
+	    !_out->about_equal (other._out, epsilon)) {
+		return false;
+	}
+
+	if (!_adjusted_white && !other._adjusted_white) {
+		return true;
+	}
+
+	if (
+		_adjusted_white && other._adjusted_white &&
+		fabs (_adjusted_white.get().x - other._adjusted_white.get().x) < epsilon &&
+		fabs (_adjusted_white.get().y - other._adjusted_white.get().y) < epsilon
+		) {
+		return true;
+	}
+
+	/* Otherwise one has an adjusted white and other hasn't, or they both have but different */
+	return false;
+}
+
+boost::numeric::ublas::matrix<double>
+ColourConversion::rgb_to_xyz () const
+{
+	/* See doc/design/colour.tex */
+
+	double const D = (_red.x - _white.x) * (_white.y - _blue.y) - (_white.x - _blue.x) * (_red.y - _white.y);
+	double const E = (_white.x - _green.x) * (_red.y - _white.y) - (_red.x - _white.x) * (_white.y - _green.y);
+	double const F = (_white.x - _green.x) * (_white.y - _blue.y) - (_white.x - _blue.x) * (_white.y - _green.y);
+	double const P = _red.y + _green.y * D / F + _blue.y * E / F;
+
+	boost::numeric::ublas::matrix<double> C (3, 3);
+	C(0, 0) = _red.x / P;
+	C(0, 1) = _green.x * D / (F * P);
+	C(0, 2) = _blue.x * E / (F * P);
+	C(1, 0) = _red.y / P;
+	C(1, 1) = _green.y * D / (F * P);
+	C(1, 2) = _blue.y * E / (F * P);
+	C(2, 0) = _red.z() / P;
+	C(2, 1) = _green.z() * D / (F * P);
+	C(2, 2) = _blue.z() * E / (F * P);
+	return C;
+}
+
+boost::numeric::ublas::matrix<double>
+ColourConversion::xyz_to_rgb () const
+{
+	boost::numeric::ublas::matrix<double> A (rgb_to_xyz ());
+
+	/* permutation matrix for the LU-factorization */
+	boost::numeric::ublas::permutation_matrix<std::size_t> pm (A.size1 ());
+
+	/* perform LU-factorization */
+	int const r = lu_factorize (A, pm);
+	DCP_ASSERT (r == 0);
+
+	/* create identity matrix of inverse */
+	boost::numeric::ublas::matrix<double> xyz_to_rgb;
+	xyz_to_rgb.assign (boost::numeric::ublas::identity_matrix<double> (A.size1 ()));
+
+	/* backsubstitute to get the inverse */
+	lu_substitute (A, pm, xyz_to_rgb);
+	
+	return xyz_to_rgb;
+}
+
+boost::numeric::ublas::matrix<double>
+ColourConversion::bradford () const
+{
+	if (!_adjusted_white || fabs (_adjusted_white.get().x) < 1e-6 || fabs (_adjusted_white.get().y) < 1e-6) {
+		boost::numeric::ublas::matrix<double> B = boost::numeric::ublas::zero_matrix<double> (3, 3);
+		B(0, 0) = 1;
+		B(1, 1) = 1;
+		B(2, 2) = 1;
+		return B;
 	}
 
-	return _in->about_equal (other._in, epsilon) && _out->about_equal (other._out, epsilon);
+	/* See doc/design/colour.tex */
+
+	boost::numeric::ublas::matrix<double> M (3, 3);
+	M(0, 0) = 0.8951;
+	M(0, 1) = 0.2664;
+	M(0, 2) = -0.1614;
+	M(1, 0) = -0.7502;
+	M(1, 1) = 1.7135;
+	M(1, 2) = 0.0367;
+	M(2, 0) = 0.0389;
+	M(2, 1) = -0.0685;
+	M(2, 2) = 1.0296;
+
+	boost::numeric::ublas::matrix<double> Mi (3, 3);
+	Mi(0, 0) = 0.9869929055;
+	Mi(0, 1) = -0.1470542564;
+	Mi(0, 2) = 0.1599626517;
+	Mi(1, 0) = 0.4323052697;
+	Mi(1, 1) = 0.5183602715;
+	Mi(1, 2) = 0.0492912282;
+	Mi(2, 0) = -0.0085286646;
+	Mi(2, 1) = 0.0400428217;
+	Mi(2, 2) = 0.9684866958;
+
+	boost::numeric::ublas::matrix<double> Gp (3, 1);
+	Gp(0, 0) = _white.x / _white.y;
+	Gp(1, 0) = 1;
+	Gp(2, 0) = (1 - _white.x - _white.y) / _white.y;
+
+	boost::numeric::ublas::matrix<double> G = boost::numeric::ublas::prod (M, Gp);
+
+	boost::numeric::ublas::matrix<double> Hp (3, 1);
+	Hp(0, 0) = _adjusted_white.get().x / _adjusted_white.get().y;
+	Hp(1, 0) = 1;
+	Hp(2, 0) = (1 - _adjusted_white.get().x - _adjusted_white.get().y) / _adjusted_white.get().y;
+
+	boost::numeric::ublas::matrix<double> H = boost::numeric::ublas::prod (M, Hp);
+
+	boost::numeric::ublas::matrix<double> C = boost::numeric::ublas::zero_matrix<double> (3, 3);
+	C(0, 0) = H(0, 0) / G(0, 0);
+	C(1, 1) = H(1, 0) / G(1, 0);
+	C(2, 2) = H(2, 0) / G(2, 0);
+
+	boost::numeric::ublas::matrix<double> CM = boost::numeric::ublas::prod (C, M);
+	return boost::numeric::ublas::prod (Mi, CM);
 }
diff --git a/src/colour_conversion.h b/src/colour_conversion.h
index cf19d447..e1ee1b11 100644
--- a/src/colour_conversion.h
+++ b/src/colour_conversion.h
@@ -1,5 +1,5 @@
 /*
-    Copyright (C) 2014 Carl Hetherington <cth@carlh.net>
+    Copyright (C) 2014-2015 Carl Hetherington <cth@carlh.net>
 
     This program is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -17,23 +17,35 @@
 
 */
 
+#include "chromaticity.h"
 #include <boost/shared_ptr.hpp>
 #include <boost/numeric/ublas/matrix.hpp>
+#include <boost/optional.hpp>
 
 namespace dcp {
 
 class TransferFunction;
 
+enum YUVToRGB {
+	YUV_TO_RGB_REC601,
+	YUV_TO_RGB_REC709,
+	YUV_TO_RGB_COUNT
+};
+
 class ColourConversion
 {
 public:
 	ColourConversion ()
-		: _matrix (3, 3)
 	{}
 	
 	ColourConversion (
 		boost::shared_ptr<const TransferFunction> in,
-		double const matrix[3][3],
+		YUVToRGB yuv_to_rgb,
+		Chromaticity red,
+		Chromaticity green,
+		Chromaticity blue,
+		Chromaticity white,
+		boost::optional<Chromaticity> adjusted_white,
 		boost::shared_ptr<const TransferFunction> out
 		);
 
@@ -41,10 +53,30 @@ public:
 		return _in;
 	}
 
-	boost::numeric::ublas::matrix<double> matrix () const {
-		return _matrix;
+	YUVToRGB yuv_to_rgb () const {
+		return _yuv_to_rgb;
+	}
+
+	Chromaticity red () const {
+		return _red;
 	}
 
+	Chromaticity green () const {
+		return _green;
+	}
+
+	Chromaticity blue () const {
+		return _blue;
+	}
+
+	Chromaticity white () const {
+		return _white;
+	}
+
+	boost::optional<Chromaticity> adjusted_white () const {
+		return _adjusted_white;
+	}
+	
 	boost::shared_ptr<const TransferFunction> out () const {
 		return _out;
 	}
@@ -53,23 +85,57 @@ public:
 		_in = f;
 	}
 
-	void set_matrix (boost::numeric::ublas::matrix<double> m) {
-		_matrix = m;
+	void set_yuv_to_rgb (YUVToRGB y) {
+		_yuv_to_rgb = y;
 	}
 
+	void set_red (Chromaticity red) {
+		_red = red;
+	}
+
+	void set_green (Chromaticity green) {
+		_green = green;
+	}
+
+	void set_blue (Chromaticity blue) {
+		_blue = blue;
+	}
+
+	void set_white (Chromaticity white) {
+		_white = white;
+	}
+
+	void set_adjusted_white (Chromaticity adjusted_white) {
+		_adjusted_white = adjusted_white;
+	}
+
+	void unset_adjusted_white () {
+		_adjusted_white = boost::optional<Chromaticity> ();
+	}
+	
 	void set_out (boost::shared_ptr<const TransferFunction> f) {
 		_out = f;
 	}
 
 	bool about_equal (ColourConversion const & other, float epsilon) const;
 
+	boost::numeric::ublas::matrix<double> rgb_to_xyz () const;
+	boost::numeric::ublas::matrix<double> xyz_to_rgb () const;
+	boost::numeric::ublas::matrix<double> bradford () const;
+
 	static ColourConversion const & srgb_to_xyz ();
-	static ColourConversion const & xyz_to_srgb ();
+	static ColourConversion const & rec601_to_xyz ();
 	static ColourConversion const & rec709_to_xyz ();
 
 protected:
 	boost::shared_ptr<const TransferFunction> _in;
-	boost::numeric::ublas::matrix<double> _matrix;
+	YUVToRGB _yuv_to_rgb;
+	Chromaticity _red;
+	Chromaticity _green;
+	Chromaticity _blue;
+	Chromaticity _white;
+	/** White point that we are adjusting to using a Bradford matrix */
+	boost::optional<Chromaticity> _adjusted_white;
 	boost::shared_ptr<const TransferFunction> _out;
 };
 
diff --git a/src/gamma_transfer_function.cc b/src/gamma_transfer_function.cc
index bf773fb3..062599d6 100644
--- a/src/gamma_transfer_function.cc
+++ b/src/gamma_transfer_function.cc
@@ -25,19 +25,18 @@ using boost::shared_ptr;
 using boost::dynamic_pointer_cast;
 using namespace dcp;
 
-GammaTransferFunction::GammaTransferFunction (bool inverse, double gamma)
-	: TransferFunction (inverse)
-	, _gamma (gamma)
+GammaTransferFunction::GammaTransferFunction (double gamma)
+	: _gamma (gamma)
 {
 
 }
 
 double *
-GammaTransferFunction::make_lut (int bit_depth) const
+GammaTransferFunction::make_lut (int bit_depth, bool inverse) const
 {
 	int const bit_length = int(std::pow(2.0f, bit_depth));
 	double* lut = new double[bit_length];
-	double const gamma = _inverse ? (1 / _gamma) : _gamma;
+	double const gamma = inverse ? (1 / _gamma) : _gamma;
 	for (int i = 0; i < bit_length; ++i) {
 		lut[i] = pow(double(i) / (bit_length - 1), gamma);
 	}
@@ -48,10 +47,6 @@ GammaTransferFunction::make_lut (int bit_depth) const
 bool
 GammaTransferFunction::about_equal (shared_ptr<const TransferFunction> other, double epsilon) const
 {
-	if (!TransferFunction::about_equal (other, epsilon)) {
-		return false;
-	}
-	
 	shared_ptr<const GammaTransferFunction> o = dynamic_pointer_cast<const GammaTransferFunction> (other);
 	if (!o) {
 		return false;
diff --git a/src/gamma_transfer_function.h b/src/gamma_transfer_function.h
index 17efe12e..282fdfe0 100644
--- a/src/gamma_transfer_function.h
+++ b/src/gamma_transfer_function.h
@@ -24,7 +24,7 @@ namespace dcp {
 class GammaTransferFunction : public TransferFunction
 {
 public:
-	GammaTransferFunction (bool inverse, double gamma);
+	GammaTransferFunction (double gamma);
 
 	double gamma () const {
 		return _gamma;
@@ -33,7 +33,7 @@ public:
 	bool about_equal (boost::shared_ptr<const TransferFunction> other, double epsilon) const;
 
 protected:
-	double * make_lut (int bit_depth) const;
+	double * make_lut (int bit_depth, bool inverse) const;
 
 private:
 	double _gamma;
diff --git a/src/modified_gamma_transfer_function.cc b/src/modified_gamma_transfer_function.cc
index fe9266bd..850c7351 100644
--- a/src/modified_gamma_transfer_function.cc
+++ b/src/modified_gamma_transfer_function.cc
@@ -25,9 +25,8 @@ using boost::shared_ptr;
 using boost::dynamic_pointer_cast;
 using namespace dcp;
 
-ModifiedGammaTransferFunction::ModifiedGammaTransferFunction (bool inverse, double power, double threshold, double A, double B)
-	: TransferFunction (inverse)
-	, _power (power)
+ModifiedGammaTransferFunction::ModifiedGammaTransferFunction (double power, double threshold, double A, double B)
+	: _power (power)
 	, _threshold (threshold)
 	, _A (A)
 	, _B (B)
@@ -36,11 +35,11 @@ ModifiedGammaTransferFunction::ModifiedGammaTransferFunction (bool inverse, doub
 }
 
 double *
-ModifiedGammaTransferFunction::make_lut (int bit_depth) const
+ModifiedGammaTransferFunction::make_lut (int bit_depth, bool inverse) const
 {
 	int const bit_length = int(std::pow(2.0f, bit_depth));
 	double* lut = new double[bit_length];
-	if (_inverse) {
+	if (inverse) {
 		double const threshold = _threshold / _B;
 		for (int i = 0; i < bit_length; ++i) {
 			double const p = static_cast<double> (i) / (bit_length - 1);
@@ -67,10 +66,6 @@ ModifiedGammaTransferFunction::make_lut (int bit_depth) const
 bool
 ModifiedGammaTransferFunction::about_equal (shared_ptr<const TransferFunction> other, double epsilon) const
 {
-	if (!TransferFunction::about_equal (other, epsilon)) {
-		return false;
-	}
-	
 	shared_ptr<const ModifiedGammaTransferFunction> o = dynamic_pointer_cast<const ModifiedGammaTransferFunction> (other);
 	if (!o) {
 		return false;
diff --git a/src/modified_gamma_transfer_function.h b/src/modified_gamma_transfer_function.h
index aec52de7..ce7303fa 100644
--- a/src/modified_gamma_transfer_function.h
+++ b/src/modified_gamma_transfer_function.h
@@ -29,7 +29,7 @@ namespace dcp {
 class ModifiedGammaTransferFunction : public TransferFunction
 {
 public:
-	ModifiedGammaTransferFunction (bool inverse, double power, double threshold, double A, double B);
+	ModifiedGammaTransferFunction (double power, double threshold, double A, double B);
 
 	double power () const {
 		return _power;
@@ -50,7 +50,7 @@ public:
 	bool about_equal (boost::shared_ptr<const TransferFunction>, double epsilon) const;
 	
 protected:
-	double * make_lut (int bit_depth) const;
+	double * make_lut (int bit_depth, bool inverse) const;
 	
 private:
 	double _power;
diff --git a/src/rgb_xyz.cc b/src/rgb_xyz.cc
index 97755b78..791d6c7c 100644
--- a/src/rgb_xyz.cc
+++ b/src/rgb_xyz.cc
@@ -70,9 +70,9 @@ dcp::xyz_to_rgba (
 	int* xyz_y = xyz_image->data (1);
 	int* xyz_z = xyz_image->data (2);
 
-	double const * lut_in = conversion.in()->lut (16);
-	double const * lut_out = conversion.out()->lut (12);
-	boost::numeric::ublas::matrix<double> matrix = conversion.matrix ();
+	double const * lut_in = conversion.in()->lut (16, true);
+	double const * lut_out = conversion.out()->lut (12, false);
+	boost::numeric::ublas::matrix<double> const matrix = conversion.xyz_to_rgb ();
 
 	int const height = xyz_image->size().height;
 	int const width = xyz_image->size().width;
@@ -150,9 +150,9 @@ dcp::xyz_to_rgb (
 	int* xyz_y = xyz_image->data (1);
 	int* xyz_z = xyz_image->data (2);
 
-	double const * lut_in = conversion.in()->lut (12);
-	double const * lut_out = conversion.out()->lut (16);
-	boost::numeric::ublas::matrix<double> matrix = conversion.matrix ();
+	double const * lut_in = conversion.in()->lut (12, true);
+	double const * lut_out = conversion.out()->lut (16, false);
+	boost::numeric::ublas::matrix<double> const matrix = conversion.xyz_to_rgb ();
 	
 	for (int y = 0; y < xyz_image->size().height; ++y) {
 		uint16_t* rgb_line = reinterpret_cast<uint16_t*> (rgb + y * stride);
@@ -238,9 +238,14 @@ dcp::rgb_to_xyz (
 		double x, y, z;
 	} d;
 
-	double const * lut_in = conversion.in()->lut (12);
-	double const * lut_out = conversion.out()->lut (16);
-	boost::numeric::ublas::matrix<double> matrix = conversion.matrix ();
+	struct {
+		double x, y, z;
+	} e;
+	
+	double const * lut_in = conversion.in()->lut (12, false);
+	double const * lut_out = conversion.out()->lut (16, true);
+	boost::numeric::ublas::matrix<double> const rgb_to_xyz = conversion.rgb_to_xyz ();
+	boost::numeric::ublas::matrix<double> const bradford = conversion.bradford ();
 
 	int jn = 0;
 	for (int y = 0; y < size.height; ++y) {
@@ -253,23 +258,27 @@ dcp::rgb_to_xyz (
 			s.b = lut_in[*p++ >> 4];
 
 			/* RGB to XYZ Matrix */
-			d.x = ((s.r * matrix(0, 0)) + (s.g * matrix(0, 1)) + (s.b * matrix(0, 2)));
-			d.y = ((s.r * matrix(1, 0)) + (s.g * matrix(1, 1)) + (s.b * matrix(1, 2)));
-			d.z = ((s.r * matrix(2, 0)) + (s.g * matrix(2, 1)) + (s.b * matrix(2, 2)));
+			d.x = ((s.r * rgb_to_xyz(0, 0)) + (s.g * rgb_to_xyz(0, 1)) + (s.b * rgb_to_xyz(0, 2)));
+			d.y = ((s.r * rgb_to_xyz(1, 0)) + (s.g * rgb_to_xyz(1, 1)) + (s.b * rgb_to_xyz(1, 2)));
+			d.z = ((s.r * rgb_to_xyz(2, 0)) + (s.g * rgb_to_xyz(2, 1)) + (s.b * rgb_to_xyz(2, 2)));
+
+			e.x = ((d.x * bradford(0, 0)) + (d.y * bradford(0, 1)) + (d.z * bradford(0, 2)));
+			e.y = ((d.x * bradford(1, 0)) + (d.y * bradford(1, 1)) + (d.z * bradford(1, 2)));
+			e.z = ((d.x * bradford(2, 0)) + (d.y * bradford(2, 1)) + (d.z * bradford(2, 2)));
 			
 			/* DCI companding */
-			d.x = d.x * DCI_COEFFICIENT * 65535;
-			d.y = d.y * DCI_COEFFICIENT * 65535;
-			d.z = d.z * DCI_COEFFICIENT * 65535;
+			e.x = e.x * DCI_COEFFICIENT * 65535;
+			e.y = e.y * DCI_COEFFICIENT * 65535;
+			e.z = e.z * DCI_COEFFICIENT * 65535;
 
-			DCP_ASSERT (d.x >= 0 && d.x < 65536);
-			DCP_ASSERT (d.y >= 0 && d.y < 65536);
-			DCP_ASSERT (d.z >= 0 && d.z < 65536);
+			DCP_ASSERT (e.x >= 0 && e.x < 65536);
+			DCP_ASSERT (e.y >= 0 && e.y < 65536);
+			DCP_ASSERT (e.z >= 0 && e.z < 65536);
 			
 			/* Out gamma LUT */
-			xyz->data(0)[jn] = lut_out[int(rint(d.x))] * 4095;
-			xyz->data(1)[jn] = lut_out[int(rint(d.y))] * 4095;
-			xyz->data(2)[jn] = lut_out[int(rint(d.z))] * 4095;
+			xyz->data(0)[jn] = lut_out[int(rint(e.x))] * 4095;
+			xyz->data(1)[jn] = lut_out[int(rint(e.y))] * 4095;
+			xyz->data(2)[jn] = lut_out[int(rint(e.z))] * 4095;
 
 			++jn;
 		}
diff --git a/src/transfer_function.cc b/src/transfer_function.cc
index 8a844d9b..bf9b0f2b 100644
--- a/src/transfer_function.cc
+++ b/src/transfer_function.cc
@@ -22,20 +22,16 @@
 
 using std::pow;
 using std::map;
+using std::pair;
+using std::make_pair;
 using boost::shared_ptr;
 using namespace dcp;
 
-TransferFunction::TransferFunction (bool inverse)
-	: _inverse (inverse)
-{
-
-}
-
 TransferFunction::~TransferFunction ()
 {
 	boost::mutex::scoped_lock lm (_mutex);
 
-	for (map<int, double*>::const_iterator i = _luts.begin(); i != _luts.end(); ++i) {
+	for (map<pair<int, bool>, double*>::const_iterator i = _luts.begin(); i != _luts.end(); ++i) {
 		delete[] i->second;
 	}
 
@@ -43,21 +39,15 @@ TransferFunction::~TransferFunction ()
 }
 
 double const *
-TransferFunction::lut (int bit_depth) const
+TransferFunction::lut (int bit_depth, bool inverse) const
 {
 	boost::mutex::scoped_lock lm (_mutex);
 	
-	map<int, double*>::const_iterator i = _luts.find (bit_depth);
+	map<pair<int, bool>, double*>::const_iterator i = _luts.find (make_pair (bit_depth, inverse));
 	if (i != _luts.end ()) {
 		return i->second;
 	}
 
-	_luts[bit_depth] = make_lut (bit_depth);
-	return _luts[bit_depth];
-}
-
-bool
-TransferFunction::about_equal (shared_ptr<const TransferFunction> other, double) const
-{
-	return _inverse == other->_inverse;
+	_luts[make_pair(bit_depth, inverse)] = make_lut (bit_depth, inverse);
+	return _luts[make_pair(bit_depth, inverse)];
 }
diff --git a/src/transfer_function.h b/src/transfer_function.h
index 4105f3d4..7f7187b5 100644
--- a/src/transfer_function.h
+++ b/src/transfer_function.h
@@ -30,22 +30,18 @@ namespace dcp {
 class TransferFunction : public boost::noncopyable
 {
 public:
-	TransferFunction (bool inverse);
-	
 	virtual ~TransferFunction ();
 
 	/** @return A look-up table (of size 2^bit_depth) whose values range from 0 to 1 */
-	double const * lut (int bit_depth) const;
+	double const * lut (int bit_depth, bool inverse) const;
 
-	virtual bool about_equal (boost::shared_ptr<const TransferFunction> other, double epsilon) const;
+	virtual bool about_equal (boost::shared_ptr<const TransferFunction> other, double epsilon) const = 0;
 
 protected:
-	virtual double * make_lut (int bit_depth) const = 0;
-
-	bool _inverse;
+	virtual double * make_lut (int bit_depth, bool inverse) const = 0;
 
 private:
-	mutable std::map<int, double*> _luts;
+	mutable std::map<std::pair<int, bool>, double*> _luts;
 	/** mutex to protect _luts */
 	mutable boost::mutex _mutex;
 };
diff --git a/src/wscript b/src/wscript
index f85b0032..1107d8cd 100644
--- a/src/wscript
+++ b/src/wscript
@@ -65,6 +65,7 @@ def build(bld):
               asset.h
               certificate_chain.h
               certificates.h
+              chromaticity.h
               colour_conversion.h
               colour_matrix.h
               cpl.h
diff --git a/test/colour_conversion_test.cc b/test/colour_conversion_test.cc
index 9111b918..2eb5edf5 100644
--- a/test/colour_conversion_test.cc
+++ b/test/colour_conversion_test.cc
@@ -28,9 +28,9 @@ using boost::shared_ptr;
 using namespace dcp;
 
 static void
-check_gamma (shared_ptr<const TransferFunction> tf, int bit_depth, float gamma)
+check_gamma (shared_ptr<const TransferFunction> tf, int bit_depth, bool inverse, float gamma)
 {
-	double const * lut = tf->lut (bit_depth);
+	double const * lut = tf->lut (bit_depth, inverse);
 	int const count = rint (pow (2.0, bit_depth));
 
 	for (int i = 0; i < count; ++i) {
@@ -39,9 +39,9 @@ check_gamma (shared_ptr<const TransferFunction> tf, int bit_depth, float gamma)
 }
 
 static void
-check_modified_gamma (shared_ptr<const TransferFunction> tf, int bit_depth, double power, double threshold, double A, double B)
+check_modified_gamma (shared_ptr<const TransferFunction> tf, int bit_depth, bool inverse, double power, double threshold, double A, double B)
 {
-	double const * lut = tf->lut (bit_depth);
+	double const * lut = tf->lut (bit_depth, inverse);
 	int const count = rint (pow (2.0, bit_depth));
 
 	for (int i = 0; i < count; ++i) {
@@ -58,25 +58,25 @@ BOOST_AUTO_TEST_CASE (colour_conversion_test1)
 {
 	ColourConversion cc = ColourConversion::srgb_to_xyz ();
 
-	check_modified_gamma (cc.in(), 8, 2.4, 0.04045, 0.055, 12.92);
-	check_modified_gamma (cc.in(), 12, 2.4, 0.04045, 0.055, 12.92);
-	check_modified_gamma (cc.in(), 16, 2.4, 0.04045, 0.055, 12.92);
+	check_modified_gamma (cc.in(), 8, false, 2.4, 0.04045, 0.055, 12.92);
+	check_modified_gamma (cc.in(), 12, false, 2.4, 0.04045, 0.055, 12.92);
+	check_modified_gamma (cc.in(), 16, false, 2.4, 0.04045, 0.055, 12.92);
 
-	check_gamma (cc.out(), 8, 1 / 2.6);
-	check_gamma (cc.out(), 12, 1 / 2.6);
-	check_gamma (cc.out(), 16, 1 / 2.6);
+	check_gamma (cc.out(), 8, true, 1 / 2.6);
+	check_gamma (cc.out(), 12, true, 1 / 2.6);
+	check_gamma (cc.out(), 16, true, 1 / 2.6);
 }
 
 BOOST_AUTO_TEST_CASE (colour_conversion_test2)
 {
 	ColourConversion cc = ColourConversion::rec709_to_xyz ();
 
-	check_modified_gamma (cc.in(), 8, 1 / 0.45, 0.081, 0.099, 4.5);
-	check_modified_gamma (cc.in(), 12, 1 / 0.45, 0.081, 0.099, 4.5);
-	check_modified_gamma (cc.in(), 16, 1 / 0.45, 0.081, 0.099, 4.5);
+	check_modified_gamma (cc.in(), 8, false, 1 / 0.45, 0.081, 0.099, 4.5);
+	check_modified_gamma (cc.in(), 12, false, 1 / 0.45, 0.081, 0.099, 4.5);
+	check_modified_gamma (cc.in(), 16, false, 1 / 0.45, 0.081, 0.099, 4.5);
 
-	check_gamma (cc.out(), 8, 1 / 2.6);
-	check_gamma (cc.out(), 12, 1 / 2.6);
-	check_gamma (cc.out(), 16, 1 / 2.6);
+	check_gamma (cc.out(), 8, true, 1 / 2.6);
+	check_gamma (cc.out(), 12, true, 1 / 2.6);
+	check_gamma (cc.out(), 16, true, 1 / 2.6);
 }
 
diff --git a/test/decryption_test.cc b/test/decryption_test.cc
index 282e74f3..0080fe3e 100644
--- a/test/decryption_test.cc
+++ b/test/decryption_test.cc
@@ -50,7 +50,7 @@ get_frame (dcp::DCP const & dcp)
 	shared_ptr<dcp::XYZImage> xyz = j2k_frame->xyz_image();
 
 	uint8_t* argb = new uint8_t[xyz->size().width * xyz->size().height * 4];
-	dcp::xyz_to_rgba (j2k_frame->xyz_image(), dcp::ColourConversion::xyz_to_srgb(), argb);
+	dcp::xyz_to_rgba (j2k_frame->xyz_image(), dcp::ColourConversion::srgb_to_xyz(), argb);
 	return make_pair (argb, xyz->size ());
 }
 
diff --git a/test/rgb_xyz_test.cc b/test/rgb_xyz_test.cc
index 26411870..0b427893 100644
--- a/test/rgb_xyz_test.cc
+++ b/test/rgb_xyz_test.cc
@@ -134,7 +134,9 @@ BOOST_AUTO_TEST_CASE (xyz_rgb_range_test)
 	scoped_array<uint8_t> rgb (new uint8_t[2 * 2 * 6]);
 
 	notes.clear ();
-	dcp::xyz_to_rgb (xyz, dcp::ColourConversion::xyz_to_srgb (), rgb.get(), 2 * 6, boost::optional<dcp::NoteHandler> (boost::bind (&note_handler, _1, _2)));
+	dcp::xyz_to_rgb (
+		xyz, dcp::ColourConversion::srgb_to_xyz (), rgb.get(), 2 * 6, boost::optional<dcp::NoteHandler> (boost::bind (&note_handler, _1, _2))
+		);
 
 	/* The 6 out-of-range samples should have been noted */
 	BOOST_REQUIRE_EQUAL (notes.size(), 6);
diff --git a/test/round_trip_test.cc b/test/round_trip_test.cc
index dcaa6939..9f1030da 100644
--- a/test/round_trip_test.cc
+++ b/test/round_trip_test.cc
@@ -110,10 +110,10 @@ BOOST_AUTO_TEST_CASE (round_trip_test)
 	shared_ptr<dcp::XYZImage> xyz_B = mxf_B->get_frame(0)->xyz_image ();
 
 	scoped_array<uint8_t> frame_A (new uint8_t[xyz_A->size().width * xyz_A->size().height * 4]);
-	dcp::xyz_to_rgba (xyz_A, dcp::ColourConversion::xyz_to_srgb(), frame_A.get());
+	dcp::xyz_to_rgba (xyz_A, dcp::ColourConversion::srgb_to_xyz(), frame_A.get());
 	
 	scoped_array<uint8_t> frame_B (new uint8_t[xyz_B->size().width * xyz_B->size().height * 4]);
-	dcp::xyz_to_rgba (xyz_B, dcp::ColourConversion::xyz_to_srgb(), frame_B.get());
+	dcp::xyz_to_rgba (xyz_B, dcp::ColourConversion::srgb_to_xyz(), frame_B.get());
 
 	BOOST_CHECK_EQUAL (xyz_A->size().width, xyz_B->size().width);
 	BOOST_CHECK_EQUAL (xyz_A->size().height, xyz_B->size().height);