Move rgb_xyz into colour_conversion.

15 files changed, 624 insertions, 772 deletions

diff --git a/benchmark/rgb_to_xyz.cc b/benchmark/rgb_to_xyz.cc
index e669e5e2..a58eb6e1 100644
--- a/benchmark/rgb_to_xyz.cc
+++ b/benchmark/rgb_to_xyz.cc
@@ -31,17 +31,20 @@
     files in the program, then also delete it here.
 */
 
-#include "openjpeg_image.h"
-#include "rgb_xyz.h"
+
 #include "colour_conversion.h"
+#include "openjpeg_image.h"
 #include <boost/scoped_array.hpp>
 #include <stdint.h>
 
+
 using boost::scoped_array;
 using std::shared_ptr;
 
+
 int const trials = 256;
 
+
 int
 main ()
 {
diff --git a/examples/read_dcp.cc b/examples/read_dcp.cc
index d88c3605..e7c8b965 100644
--- a/examples/read_dcp.cc
+++ b/examples/read_dcp.cc
@@ -35,7 +35,6 @@
 #include "text_asset.h"
 #include "openjpeg_image.h"
 #include "colour_conversion.h"
-#include "rgb_xyz.h"
 /* This DISABLE/ENABLE pair is just to ignore some warnings from Magick++.h; they
  * can be removed.
  */
diff --git a/src/colour_conversion.cc b/src/colour_conversion.cc
index 3b98ccd0..8a3c8bff 100644
--- a/src/colour_conversion.cc
+++ b/src/colour_conversion.cc
@@ -38,22 +38,29 @@
 
 
 #include "colour_conversion.h"
+#include "compose.hpp"
+#include "dcp_assert.h"
 #include "gamma_transfer_function.h"
+#include "identity_transfer_function.h"
 #include "modified_gamma_transfer_function.h"
+#include "openjpeg_image.h"
 #include "s_gamut3_transfer_function.h"
-#include "identity_transfer_function.h"
-#include "dcp_assert.h"
 #include <boost/numeric/ublas/matrix.hpp>
 #include <boost/numeric/ublas/lu.hpp>
 #include <boost/numeric/ublas/io.hpp>
 
 
-using std::shared_ptr;
 using std::make_shared;
+using std::max;
+using std::min;
+using std::shared_ptr;
 using boost::optional;
 using namespace dcp;
 
 
+static auto constexpr DCI_COEFFICIENT = 48.0 / 52.37;
+
+
 ColourConversion const &
 ColourConversion::srgb_to_dcp()
 {
@@ -343,3 +350,294 @@ ColourConversion::bradford() const
 	boost::numeric::ublas::matrix<double> CM = boost::numeric::ublas::prod(C, M);
 	return boost::numeric::ublas::prod(Mi, CM);
 }
+
+
+
+
+void
+dcp::xyz_to_rgba(
+	std::shared_ptr<const OpenJPEGImage> xyz_image,
+	ColourConversion const & conversion,
+	uint8_t* argb,
+	int stride
+	)
+{
+	int const max_colour = pow(2, 16) - 1;
+
+	struct {
+		double x, y, z;
+	} s;
+
+	struct {
+		double r, g, b;
+	} d;
+
+	int* xyz_x = xyz_image->data(0);
+	int* xyz_y = xyz_image->data(1);
+	int* xyz_z = xyz_image->data(2);
+
+	auto lut_in = conversion.out_j2k()->double_lut(0, 1, 12, false);
+	auto lut_out = conversion.in()->double_lut(0, 1, 16, true);
+	boost::numeric::ublas::matrix<double> const matrix = conversion.xyz_to_rgb();
+
+	double fast_matrix[9] = {
+		matrix(0, 0), matrix(0, 1), matrix(0, 2),
+		matrix(1, 0), matrix(1, 1), matrix(1, 2),
+		matrix(2, 0), matrix(2, 1), matrix(2, 2)
+	};
+
+	int const height = xyz_image->size().height;
+	int const width = xyz_image->size().width;
+
+	for (int y = 0; y < height; ++y) {
+		uint8_t* argb_line = argb;
+		for (int x = 0; x < width; ++x) {
+
+			DCP_ASSERT(*xyz_x >= 0 && *xyz_y >= 0 && *xyz_z >= 0 && *xyz_x < 4096 && *xyz_y < 4096 && *xyz_z < 4096);
+
+			/* In gamma LUT */
+			s.x = lut_in[*xyz_x++];
+			s.y = lut_in[*xyz_y++];
+			s.z = lut_in[*xyz_z++];
+
+			/* DCI companding */
+			s.x /= DCI_COEFFICIENT;
+			s.y /= DCI_COEFFICIENT;
+			s.z /= DCI_COEFFICIENT;
+
+			/* XYZ to RGB */
+			d.r = ((s.x * fast_matrix[0]) + (s.y * fast_matrix[1]) + (s.z * fast_matrix[2]));
+			d.g = ((s.x * fast_matrix[3]) + (s.y * fast_matrix[4]) + (s.z * fast_matrix[5]));
+			d.b = ((s.x * fast_matrix[6]) + (s.y * fast_matrix[7]) + (s.z * fast_matrix[8]));
+
+			d.r = min(d.r, 1.0);
+			d.r = max(d.r, 0.0);
+
+			d.g = min(d.g, 1.0);
+			d.g = max(d.g, 0.0);
+
+			d.b = min(d.b, 1.0);
+			d.b = max(d.b, 0.0);
+
+			/* Out gamma LUT */
+			*argb_line++ = lut_out[lrint(d.b * max_colour)] * 0xff;
+			*argb_line++ = lut_out[lrint(d.g * max_colour)] * 0xff;
+			*argb_line++ = lut_out[lrint(d.r * max_colour)] * 0xff;
+			*argb_line++ = 0xff;
+		}
+
+		argb += stride;
+	}
+}
+
+
+void
+dcp::xyz_to_rgb(
+	shared_ptr<const OpenJPEGImage> xyz_image,
+	ColourConversion const & conversion,
+	uint8_t* rgb,
+	int stride,
+	optional<NoteHandler> note
+	)
+{
+	struct {
+		double x, y, z;
+	} s;
+
+	struct {
+		double r, g, b;
+	} d;
+
+	/* These should be 12-bit values from 0-4095 */
+	int* xyz_x = xyz_image->data(0);
+	int* xyz_y = xyz_image->data(1);
+	int* xyz_z = xyz_image->data(2);
+
+	auto lut_in = conversion.out_j2k()->double_lut(0, 1, 12, false);
+	auto lut_out = conversion.in()->double_lut(0, 1, 16, true);
+	auto const matrix = conversion.xyz_to_rgb();
+
+	double fast_matrix[9] = {
+		matrix(0, 0), matrix(0, 1), matrix(0, 2),
+		matrix(1, 0), matrix(1, 1), matrix(1, 2),
+		matrix(2, 0), matrix(2, 1), matrix(2, 2)
+	};
+
+	int const height = xyz_image->size().height;
+	int const width = xyz_image->size().width;
+
+	for (int y = 0; y < height; ++y) {
+		auto rgb_line = reinterpret_cast<uint16_t*>(rgb + y * stride);
+		for (int x = 0; x < width; ++x) {
+
+			int cx = *xyz_x++;
+			int cy = *xyz_y++;
+			int cz = *xyz_z++;
+
+			if (cx < 0 || cx > 4095) {
+				if (note) {
+					note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cx));
+				}
+				cx = max(min(cx, 4095), 0);
+			}
+
+			if (cy < 0 || cy > 4095) {
+				if (note) {
+					note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cy));
+				}
+				cy = max(min(cy, 4095), 0);
+			}
+
+			if (cz < 0 || cz > 4095) {
+				if (note) {
+					note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cz));
+				}
+				cz = max(min(cz, 4095), 0);
+			}
+
+			/* In gamma LUT */
+			s.x = lut_in[cx];
+			s.y = lut_in[cy];
+			s.z = lut_in[cz];
+
+			/* DCI companding */
+			s.x /= DCI_COEFFICIENT;
+			s.y /= DCI_COEFFICIENT;
+			s.z /= DCI_COEFFICIENT;
+
+			/* XYZ to RGB */
+			d.r = ((s.x * fast_matrix[0]) + (s.y * fast_matrix[1]) + (s.z * fast_matrix[2]));
+			d.g = ((s.x * fast_matrix[3]) + (s.y * fast_matrix[4]) + (s.z * fast_matrix[5]));
+			d.b = ((s.x * fast_matrix[6]) + (s.y * fast_matrix[7]) + (s.z * fast_matrix[8]));
+
+			d.r = min(d.r, 1.0);
+			d.r = max(d.r, 0.0);
+
+			d.g = min(d.g, 1.0);
+			d.g = max(d.g, 0.0);
+
+			d.b = min(d.b, 1.0);
+			d.b = max(d.b, 0.0);
+
+			*rgb_line++ = lrint(lut_out[lrint(d.r * 65535)] * 65535);
+			*rgb_line++ = lrint(lut_out[lrint(d.g * 65535)] * 65535);
+			*rgb_line++ = lrint(lut_out[lrint(d.b * 65535)] * 65535);
+		}
+	}
+}
+
+void
+dcp::combined_rgb_to_xyz(ColourConversion const & conversion, double* matrix)
+{
+	auto const rgb_to_xyz = conversion.rgb_to_xyz();
+	auto const bradford = conversion.bradford();
+
+	matrix[0] = (bradford(0, 0) * rgb_to_xyz(0, 0) + bradford(0, 1) * rgb_to_xyz(1, 0) + bradford(0, 2) * rgb_to_xyz(2, 0)) * DCI_COEFFICIENT;
+	matrix[1] = (bradford(0, 0) * rgb_to_xyz(0, 1) + bradford(0, 1) * rgb_to_xyz(1, 1) + bradford(0, 2) * rgb_to_xyz(2, 1)) * DCI_COEFFICIENT;
+	matrix[2] = (bradford(0, 0) * rgb_to_xyz(0, 2) + bradford(0, 1) * rgb_to_xyz(1, 2) + bradford(0, 2) * rgb_to_xyz(2, 2)) * DCI_COEFFICIENT;
+	matrix[3] = (bradford(1, 0) * rgb_to_xyz(0, 0) + bradford(1, 1) * rgb_to_xyz(1, 0) + bradford(1, 2) * rgb_to_xyz(2, 0)) * DCI_COEFFICIENT;
+	matrix[4] = (bradford(1, 0) * rgb_to_xyz(0, 1) + bradford(1, 1) * rgb_to_xyz(1, 1) + bradford(1, 2) * rgb_to_xyz(2, 1)) * DCI_COEFFICIENT;
+	matrix[5] = (bradford(1, 0) * rgb_to_xyz(0, 2) + bradford(1, 1) * rgb_to_xyz(1, 2) + bradford(1, 2) * rgb_to_xyz(2, 2)) * DCI_COEFFICIENT;
+	matrix[6] = (bradford(2, 0) * rgb_to_xyz(0, 0) + bradford(2, 1) * rgb_to_xyz(1, 0) + bradford(2, 2) * rgb_to_xyz(2, 0)) * DCI_COEFFICIENT;
+	matrix[7] = (bradford(2, 0) * rgb_to_xyz(0, 1) + bradford(2, 1) * rgb_to_xyz(1, 1) + bradford(2, 2) * rgb_to_xyz(2, 1)) * DCI_COEFFICIENT;
+	matrix[8] = (bradford(2, 0) * rgb_to_xyz(0, 2) + bradford(2, 1) * rgb_to_xyz(1, 2) + bradford(2, 2) * rgb_to_xyz(2, 2)) * DCI_COEFFICIENT;
+}
+
+
+PiecewiseLUT2
+dcp::make_inverse_gamma_lut(shared_ptr<const TransferFunction> fn)
+{
+	/* The parameters here were chosen by trial and error to reduce errors when running rgb_xyz_lut_test */
+	return PiecewiseLUT2(fn, 0.062, 16, 12, true, 4095);
+}
+
+
+template <class T>
+void
+rgb_to_xyz_internal(
+	uint8_t const* rgb,
+	T*& xyz_x,
+	T*& xyz_y,
+	T*& xyz_z,
+	dcp::Size size,
+	int stride,
+	ColourConversion const& conversion
+	)
+{
+	struct {
+		double r, g, b;
+	} s;
+
+	struct {
+		double x, y, z;
+	} d;
+
+	auto lut_in = conversion.in()->double_lut(0, 1, 12, false);
+	auto lut_out = make_inverse_gamma_lut(conversion.out_j2k());
+
+	/* This is is the product of the RGB to XYZ matrix, the Bradford transform and the DCI companding */
+	double fast_matrix[9];
+	combined_rgb_to_xyz(conversion, fast_matrix);
+
+	for (int y = 0; y < size.height; ++y) {
+		auto p = reinterpret_cast<uint16_t const *>(rgb + y * stride);
+		for (int x = 0; x < size.width; ++x) {
+
+			/* In gamma LUT (converting 16-bit to 12-bit) */
+			s.r = lut_in[*p++ >> 4];
+			s.g = lut_in[*p++ >> 4];
+			s.b = lut_in[*p++ >> 4];
+
+			/* RGB to XYZ, Bradford transform and DCI companding */
+			d.x = s.r * fast_matrix[0] + s.g * fast_matrix[1] + s.b * fast_matrix[2];
+			d.y = s.r * fast_matrix[3] + s.g * fast_matrix[4] + s.b * fast_matrix[5];
+			d.z = s.r * fast_matrix[6] + s.g * fast_matrix[7] + s.b * fast_matrix[8];
+
+			/* Clamp */
+			d.x = max(0.0, d.x);
+			d.y = max(0.0, d.y);
+			d.z = max(0.0, d.z);
+			d.x = min(1.0, d.x);
+			d.y = min(1.0, d.y);
+			d.z = min(1.0, d.z);
+
+			/* Out gamma LUT */
+			*xyz_x++ = lut_out.lookup(d.x);
+			*xyz_y++ = lut_out.lookup(d.y);
+			*xyz_z++ = lut_out.lookup(d.z);
+		}
+	}
+}
+
+
+shared_ptr<dcp::OpenJPEGImage>
+dcp::rgb_to_xyz(
+	uint8_t const * rgb,
+	dcp::Size size,
+	int stride,
+	ColourConversion const & conversion
+	)
+{
+	auto xyz = make_shared<OpenJPEGImage>(size);
+
+	int* xyz_x = xyz->data(0);
+	int* xyz_y = xyz->data(1);
+	int* xyz_z = xyz->data(2);
+
+	rgb_to_xyz_internal(rgb, xyz_x, xyz_y, xyz_z, size, stride, conversion);
+
+	return xyz;
+}
+
+
+void
+dcp::rgb_to_xyz(
+	uint8_t const * rgb,
+	uint16_t* dst,
+	dcp::Size size,
+	int stride,
+	ColourConversion const & conversion
+	)
+{
+	rgb_to_xyz_internal(rgb, dst, dst, dst, size, stride, conversion);
+}
diff --git a/src/colour_conversion.h b/src/colour_conversion.h
index d96a66b9..3c39b45a 100644
--- a/src/colour_conversion.h
+++ b/src/colour_conversion.h
@@ -42,17 +42,22 @@
 
 
 #include "chromaticity.h"
+#include "piecewise_lut.h"
+#include "types.h"
 #include <memory>
 #if BOOST_VERSION >= 106400
 #include <boost/serialization/array_wrapper.hpp>
 #endif
 #include <boost/numeric/ublas/matrix.hpp>
 #include <boost/optional.hpp>
+#include <stdint.h>
 
 
 namespace dcp {
 
 
+class Image;
+class OpenJPEGImage;
 class TransferFunction;
 
 
@@ -199,7 +204,87 @@ protected:
 };
 
 
-}
+/** Convert an XYZ image to RGBA.
+ *  @param xyz_image Image in XYZ.
+ *  @param conversion Colour conversion to use.
+ *  @param argb Buffer to fill with RGBA data.  The format of the data is:
+ *
+ *  <pre>
+ *  Byte   /- 0 -------|- 1 --------|- 2 --------|- 3 --------|- 4 --------|- 5 --------| ...
+ *         |(0, 0) Blue|(0, 0)Green |(0, 0) Red  |(0, 0) Alpha|(0, 1) Blue |(0, 1) Green| ...
+ *  </pre>
+ *
+ *  So that the first byte is the blue component of the pixel at x=0, y=0, the second
+ *  is the green component, and so on.
+ *
+ *  Lines are packed so that the second row directly follows the first.
+ */
+extern void xyz_to_rgba(
+	std::shared_ptr<const OpenJPEGImage>,
+	ColourConversion const & conversion,
+	uint8_t* rgba,
+	int stride
+	);
+
+
+/** Convert an XYZ image to 48bpp RGB.
+ *  @param xyz_image Frame in XYZ.
+ *  @param conversion Colour conversion to use.
+ *  @param rgb Buffer to fill with RGB data.  Format is packed RGB
+ *  16:16:16, 48bpp, 16R, 16G, 16B, with the 2-byte value for each
+ *  R/G/B component stored as little-endian; i.e. AV_PIX_FMT_RGB48LE.
+ *  @param stride Stride for RGB data in bytes.
+ *  @param note Optional handler for any notes that may be made during the conversion (e.g. when clamping occurs).
+ */
+extern void xyz_to_rgb(
+	std::shared_ptr<const OpenJPEGImage>,
+	ColourConversion const & conversion,
+	uint8_t* rgb,
+	int stride,
+	boost::optional<NoteHandler> note = boost::optional<NoteHandler> ()
+	);
+
+
+extern PiecewiseLUT2 make_inverse_gamma_lut(std::shared_ptr<const TransferFunction> fn);
+
+
+/** @param rgb RGB data; packed RGB 16:16:16, 48bpp, 16R, 16G, 16B,
+ *  with the 2-byte value for each R/G/B component stored as
+ *  little-endian; i.e. AV_PIX_FMT_RGB48LE.
+ *  @param dst Buffer to fill with packed 16-bit XYZ data, i.e. first
+ *  16-bit word is X, second is Y etc.
+ *  @param size size of RGB image in pixels.
+ *  @param size stride of RGB data in pixels.
+ */
+extern void rgb_to_xyz(
+	uint8_t const * rgb,
+	uint16_t* dst,
+	dcp::Size size,
+	int stride,
+	ColourConversion const& conversion
+	);
+
+
+/** @param rgb RGB data; packed RGB 16:16:16, 48bpp, 16R, 16G, 16B,
+ *  with the 2-byte value for each R/G/B component stored as
+ *  little-endian; i.e. AV_PIX_FMT_RGB48LE.
+ *  @param size size of RGB image in pixels.
+ *  @param size stride of RGB data in pixels.
+ */
+extern std::shared_ptr<OpenJPEGImage> rgb_to_xyz(
+	uint8_t const * rgb,
+	dcp::Size size,
+	int stride,
+	ColourConversion const& conversion
+	);
 
 
+/** @param conversion Colour conversion.
+ *  @param matrix Filled in with the product of the RGB to XYZ matrix, the Bradford transform and the DCI companding.
+ */
+extern void combined_rgb_to_xyz(ColourConversion const & conversion, double* matrix);
+
+
+}
+
 #endif
diff --git a/src/mono_j2k_picture_frame.cc b/src/mono_j2k_picture_frame.cc
index 69ba9a59..8bc07a6d 100644
--- a/src/mono_j2k_picture_frame.cc
+++ b/src/mono_j2k_picture_frame.cc
@@ -45,7 +45,6 @@
 #include "filesystem.h"
 #include "j2k_transcode.h"
 #include "mono_j2k_picture_frame.h"
-#include "rgb_xyz.h"
 #include "util.h"
 #include <asdcp/KM_fileio.h>
 #include <asdcp/AS_DCP.h>
diff --git a/src/rgb_xyz.cc b/src/rgb_xyz.cc
deleted file mode 100644
index eaf1931b..00000000
--- a/src/rgb_xyz.cc
+++ /dev/null
@@ -1,357 +0,0 @@
-/*
-    Copyright (C) 2013-2021 Carl Hetherington <cth@carlh.net>
-
-    This file is part of libdcp.
-
-    libdcp 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.
-
-    libdcp 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 libdcp.  If not, see <http://www.gnu.org/licenses/>.
-
-    In addition, as a special exception, the copyright holders give
-    permission to link the code of portions of this program with the
-    OpenSSL library under certain conditions as described in each
-    individual source file, and distribute linked combinations
-    including the two.
-
-    You must obey the GNU General Public License in all respects
-    for all of the code used other than OpenSSL.  If you modify
-    file(s) with this exception, you may extend this exception to your
-    version of the file(s), but you are not obligated to do so.  If you
-    do not wish to do so, delete this exception statement from your
-    version.  If you delete this exception statement from all source
-    files in the program, then also delete it here.
-*/
-
-
-/** @file  rgb_xyz.cc
- *  @brief Conversion between RGB and XYZ
- */
-
-
-#include "colour_conversion.h"
-#include "compose.hpp"
-#include "dcp_assert.h"
-#include "openjpeg_image.h"
-#include "piecewise_lut.h"
-#include "rgb_xyz.h"
-#include "transfer_function.h"
-#include <cmath>
-
-
-using std::cout;
-using std::make_shared;
-using std::max;
-using std::min;
-using std::shared_ptr;
-using boost::optional;
-using namespace dcp;
-
-
-static auto constexpr DCI_COEFFICIENT = 48.0 / 52.37;
-
-
-void
-dcp::xyz_to_rgba (
-	std::shared_ptr<const OpenJPEGImage> xyz_image,
-	ColourConversion const & conversion,
-	uint8_t* argb,
-	int stride
-	)
-{
-	int const max_colour = pow (2, 16) - 1;
-
-	struct {
-		double x, y, z;
-	} s;
-
-	struct {
-		double r, g, b;
-	} d;
-
-	int* xyz_x = xyz_image->data (0);
-	int* xyz_y = xyz_image->data (1);
-	int* xyz_z = xyz_image->data (2);
-
-	auto lut_in = conversion.out_j2k()->double_lut(0, 1, 12, false);
-	auto lut_out = conversion.in()->double_lut(0, 1, 16, true);
-	boost::numeric::ublas::matrix<double> const matrix = conversion.xyz_to_rgb ();
-
-	double fast_matrix[9] = {
-		matrix (0, 0), matrix (0, 1), matrix (0, 2),
-		matrix (1, 0), matrix (1, 1), matrix (1, 2),
-		matrix (2, 0), matrix (2, 1), matrix (2, 2)
-	};
-
-	int const height = xyz_image->size().height;
-	int const width = xyz_image->size().width;
-
-	for (int y = 0; y < height; ++y) {
-		uint8_t* argb_line = argb;
-		for (int x = 0; x < width; ++x) {
-
-			DCP_ASSERT (*xyz_x >= 0 && *xyz_y >= 0 && *xyz_z >= 0 && *xyz_x < 4096 && *xyz_y < 4096 && *xyz_z < 4096);
-
-			/* In gamma LUT */
-			s.x = lut_in[*xyz_x++];
-			s.y = lut_in[*xyz_y++];
-			s.z = lut_in[*xyz_z++];
-
-			/* DCI companding */
-			s.x /= DCI_COEFFICIENT;
-			s.y /= DCI_COEFFICIENT;
-			s.z /= DCI_COEFFICIENT;
-
-			/* XYZ to RGB */
-			d.r = ((s.x * fast_matrix[0]) + (s.y * fast_matrix[1]) + (s.z * fast_matrix[2]));
-			d.g = ((s.x * fast_matrix[3]) + (s.y * fast_matrix[4]) + (s.z * fast_matrix[5]));
-			d.b = ((s.x * fast_matrix[6]) + (s.y * fast_matrix[7]) + (s.z * fast_matrix[8]));
-
-			d.r = min (d.r, 1.0);
-			d.r = max (d.r, 0.0);
-
-			d.g = min (d.g, 1.0);
-			d.g = max (d.g, 0.0);
-
-			d.b = min (d.b, 1.0);
-			d.b = max (d.b, 0.0);
-
-			/* Out gamma LUT */
-			*argb_line++ = lut_out[lrint(d.b * max_colour)] * 0xff;
-			*argb_line++ = lut_out[lrint(d.g * max_colour)] * 0xff;
-			*argb_line++ = lut_out[lrint(d.r * max_colour)] * 0xff;
-			*argb_line++ = 0xff;
-		}
-
-		argb += stride;
-	}
-}
-
-
-void
-dcp::xyz_to_rgb (
-	shared_ptr<const OpenJPEGImage> xyz_image,
-	ColourConversion const & conversion,
-	uint8_t* rgb,
-	int stride,
-	optional<NoteHandler> note
-	)
-{
-	struct {
-		double x, y, z;
-	} s;
-
-	struct {
-		double r, g, b;
-	} d;
-
-	/* These should be 12-bit values from 0-4095 */
-	int* xyz_x = xyz_image->data (0);
-	int* xyz_y = xyz_image->data (1);
-	int* xyz_z = xyz_image->data (2);
-
-	auto lut_in = conversion.out_j2k()->double_lut(0, 1, 12, false);
-	auto lut_out = conversion.in()->double_lut(0, 1, 16, true);
-	auto const matrix = conversion.xyz_to_rgb ();
-
-	double fast_matrix[9] = {
-		matrix (0, 0), matrix (0, 1), matrix (0, 2),
-		matrix (1, 0), matrix (1, 1), matrix (1, 2),
-		matrix (2, 0), matrix (2, 1), matrix (2, 2)
-	};
-
-	int const height = xyz_image->size().height;
-	int const width = xyz_image->size().width;
-
-	for (int y = 0; y < height; ++y) {
-		auto rgb_line = reinterpret_cast<uint16_t*> (rgb + y * stride);
-		for (int x = 0; x < width; ++x) {
-
-			int cx = *xyz_x++;
-			int cy = *xyz_y++;
-			int cz = *xyz_z++;
-
-			if (cx < 0 || cx > 4095) {
-				if (note) {
-					note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cx));
-				}
-				cx = max (min (cx, 4095), 0);
-			}
-
-			if (cy < 0 || cy > 4095) {
-				if (note) {
-					note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cy));
-				}
-				cy = max (min (cy, 4095), 0);
-			}
-
-			if (cz < 0 || cz > 4095) {
-				if (note) {
-					note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cz));
-				}
-				cz = max (min (cz, 4095), 0);
-			}
-
-			/* In gamma LUT */
-			s.x = lut_in[cx];
-			s.y = lut_in[cy];
-			s.z = lut_in[cz];
-
-			/* DCI companding */
-			s.x /= DCI_COEFFICIENT;
-			s.y /= DCI_COEFFICIENT;
-			s.z /= DCI_COEFFICIENT;
-
-			/* XYZ to RGB */
-			d.r = ((s.x * fast_matrix[0]) + (s.y * fast_matrix[1]) + (s.z * fast_matrix[2]));
-			d.g = ((s.x * fast_matrix[3]) + (s.y * fast_matrix[4]) + (s.z * fast_matrix[5]));
-			d.b = ((s.x * fast_matrix[6]) + (s.y * fast_matrix[7]) + (s.z * fast_matrix[8]));
-
-			d.r = min (d.r, 1.0);
-			d.r = max (d.r, 0.0);
-
-			d.g = min (d.g, 1.0);
-			d.g = max (d.g, 0.0);
-
-			d.b = min (d.b, 1.0);
-			d.b = max (d.b, 0.0);
-
-			*rgb_line++ = lrint(lut_out[lrint(d.r * 65535)] * 65535);
-			*rgb_line++ = lrint(lut_out[lrint(d.g * 65535)] * 65535);
-			*rgb_line++ = lrint(lut_out[lrint(d.b * 65535)] * 65535);
-		}
-	}
-}
-
-void
-dcp::combined_rgb_to_xyz (ColourConversion const & conversion, double* matrix)
-{
-	auto const rgb_to_xyz = conversion.rgb_to_xyz ();
-	auto const bradford = conversion.bradford ();
-
-	matrix[0] = (bradford (0, 0) * rgb_to_xyz (0, 0) + bradford (0, 1) * rgb_to_xyz (1, 0) + bradford (0, 2) * rgb_to_xyz (2, 0))
-		* DCI_COEFFICIENT;
-	matrix[1] = (bradford (0, 0) * rgb_to_xyz (0, 1) + bradford (0, 1) * rgb_to_xyz (1, 1) + bradford (0, 2) * rgb_to_xyz (2, 1))
-		* DCI_COEFFICIENT;
-	matrix[2] = (bradford (0, 0) * rgb_to_xyz (0, 2) + bradford (0, 1) * rgb_to_xyz (1, 2) + bradford (0, 2) * rgb_to_xyz (2, 2))
-		* DCI_COEFFICIENT;
-	matrix[3] = (bradford (1, 0) * rgb_to_xyz (0, 0) + bradford (1, 1) * rgb_to_xyz (1, 0) + bradford (1, 2) * rgb_to_xyz (2, 0))
-		* DCI_COEFFICIENT;
-	matrix[4] = (bradford (1, 0) * rgb_to_xyz (0, 1) + bradford (1, 1) * rgb_to_xyz (1, 1) + bradford (1, 2) * rgb_to_xyz (2, 1))
-		* DCI_COEFFICIENT;
-	matrix[5] = (bradford (1, 0) * rgb_to_xyz (0, 2) + bradford (1, 1) * rgb_to_xyz (1, 2) + bradford (1, 2) * rgb_to_xyz (2, 2))
-		* DCI_COEFFICIENT;
-	matrix[6] = (bradford (2, 0) * rgb_to_xyz (0, 0) + bradford (2, 1) * rgb_to_xyz (1, 0) + bradford (2, 2) * rgb_to_xyz (2, 0))
-		* DCI_COEFFICIENT;
-	matrix[7] = (bradford (2, 0) * rgb_to_xyz (0, 1) + bradford (2, 1) * rgb_to_xyz (1, 1) + bradford (2, 2) * rgb_to_xyz (2, 1))
-		* DCI_COEFFICIENT;
-	matrix[8] = (bradford (2, 0) * rgb_to_xyz (0, 2) + bradford (2, 1) * rgb_to_xyz (1, 2) + bradford (2, 2) * rgb_to_xyz (2, 2))
-		* DCI_COEFFICIENT;
-}
-
-
-PiecewiseLUT2
-dcp::make_inverse_gamma_lut(shared_ptr<const TransferFunction> fn)
-{
-	/* The parameters here were chosen by trial and error to reduce errors when running rgb_xyz_lut_test */
-	return PiecewiseLUT2(fn, 0.062, 16, 12, true, 4095);
-}
-
-
-template <class T>
-void
-rgb_to_xyz_internal(
-	uint8_t const* rgb,
-	T*& xyz_x,
-	T*& xyz_y,
-	T*& xyz_z,
-	dcp::Size size,
-	int stride,
-	ColourConversion const& conversion
-	)
-{
-	struct {
-		double r, g, b;
-	} s;
-
-	struct {
-		double x, y, z;
-	} d;
-
-	auto lut_in = conversion.in()->double_lut(0, 1, 12, false);
-	auto lut_out = make_inverse_gamma_lut(conversion.out_j2k());
-
-	/* This is is the product of the RGB to XYZ matrix, the Bradford transform and the DCI companding */
-	double fast_matrix[9];
-	combined_rgb_to_xyz (conversion, fast_matrix);
-
-	for (int y = 0; y < size.height; ++y) {
-		auto p = reinterpret_cast<uint16_t const *> (rgb + y * stride);
-		for (int x = 0; x < size.width; ++x) {
-
-			/* In gamma LUT (converting 16-bit to 12-bit) */
-			s.r = lut_in[*p++ >> 4];
-			s.g = lut_in[*p++ >> 4];
-			s.b = lut_in[*p++ >> 4];
-
-			/* RGB to XYZ, Bradford transform and DCI companding */
-			d.x = s.r * fast_matrix[0] + s.g * fast_matrix[1] + s.b * fast_matrix[2];
-			d.y = s.r * fast_matrix[3] + s.g * fast_matrix[4] + s.b * fast_matrix[5];
-			d.z = s.r * fast_matrix[6] + s.g * fast_matrix[7] + s.b * fast_matrix[8];
-
-			/* Clamp */
-			d.x = max (0.0, d.x);
-			d.y = max (0.0, d.y);
-			d.z = max (0.0, d.z);
-			d.x = min (1.0, d.x);
-			d.y = min (1.0, d.y);
-			d.z = min (1.0, d.z);
-
-			/* Out gamma LUT */
-			*xyz_x++ = lut_out.lookup(d.x);
-			*xyz_y++ = lut_out.lookup(d.y);
-			*xyz_z++ = lut_out.lookup(d.z);
-		}
-	}
-}
-
-
-shared_ptr<dcp::OpenJPEGImage>
-dcp::rgb_to_xyz (
-	uint8_t const * rgb,
-	dcp::Size size,
-	int stride,
-	ColourConversion const & conversion
-	)
-{
-	auto xyz = make_shared<OpenJPEGImage>(size);
-
-	int* xyz_x = xyz->data (0);
-	int* xyz_y = xyz->data (1);
-	int* xyz_z = xyz->data (2);
-
-	rgb_to_xyz_internal(rgb, xyz_x, xyz_y, xyz_z, size, stride, conversion);
-
-	return xyz;
-}
-
-
-void
-dcp::rgb_to_xyz (
-	uint8_t const * rgb,
-	uint16_t* dst,
-	dcp::Size size,
-	int stride,
-	ColourConversion const & conversion
-	)
-{
-	rgb_to_xyz_internal(rgb, dst, dst, dst, size, stride, conversion);
-}
diff --git a/src/rgb_xyz.h b/src/rgb_xyz.h
deleted file mode 100644
index 315295ea..00000000
--- a/src/rgb_xyz.h
+++ /dev/null
@@ -1,132 +0,0 @@
-/*
-    Copyright (C) 2013-2021 Carl Hetherington <cth@carlh.net>
-
-    This file is part of libdcp.
-
-    libdcp 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.
-
-    libdcp 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 libdcp.  If not, see <http://www.gnu.org/licenses/>.
-
-    In addition, as a special exception, the copyright holders give
-    permission to link the code of portions of this program with the
-    OpenSSL library under certain conditions as described in each
-    individual source file, and distribute linked combinations
-    including the two.
-
-    You must obey the GNU General Public License in all respects
-    for all of the code used other than OpenSSL.  If you modify
-    file(s) with this exception, you may extend this exception to your
-    version of the file(s), but you are not obligated to do so.  If you
-    do not wish to do so, delete this exception statement from your
-    version.  If you delete this exception statement from all source
-    files in the program, then also delete it here.
-*/
-
-
-/** @file  src/rgb_xyz.h
- *  @brief Conversion between RGB and XYZ
- */
-
-
-#include "piecewise_lut.h"
-#include <memory>
-#include <boost/optional.hpp>
-#include <stdint.h>
-
-
-namespace dcp {
-
-
-class OpenJPEGImage;
-class Image;
-class ColourConversion;
-
-
-/** Convert an XYZ image to RGBA.
- *  @param xyz_image Image in XYZ.
- *  @param conversion Colour conversion to use.
- *  @param argb Buffer to fill with RGBA data.  The format of the data is:
- *
- *  <pre>
- *  Byte   /- 0 -------|- 1 --------|- 2 --------|- 3 --------|- 4 --------|- 5 --------| ...
- *         |(0, 0) Blue|(0, 0)Green |(0, 0) Red  |(0, 0) Alpha|(0, 1) Blue |(0, 1) Green| ...
- *  </pre>
- *
- *  So that the first byte is the blue component of the pixel at x=0, y=0, the second
- *  is the green component, and so on.
- *
- *  Lines are packed so that the second row directly follows the first.
- */
-extern void xyz_to_rgba (
-	std::shared_ptr<const OpenJPEGImage>,
-	ColourConversion const & conversion,
-	uint8_t* rgba,
-	int stride
-	);
-
-
-/** Convert an XYZ image to 48bpp RGB.
- *  @param xyz_image Frame in XYZ.
- *  @param conversion Colour conversion to use.
- *  @param rgb Buffer to fill with RGB data.  Format is packed RGB
- *  16:16:16, 48bpp, 16R, 16G, 16B, with the 2-byte value for each
- *  R/G/B component stored as little-endian; i.e. AV_PIX_FMT_RGB48LE.
- *  @param stride Stride for RGB data in bytes.
- *  @param note Optional handler for any notes that may be made during the conversion (e.g. when clamping occurs).
- */
-extern void xyz_to_rgb (
-	std::shared_ptr<const OpenJPEGImage>,
-	ColourConversion const & conversion,
-	uint8_t* rgb,
-	int stride,
-	boost::optional<NoteHandler> note = boost::optional<NoteHandler> ()
-	);
-
-
-extern PiecewiseLUT2 make_inverse_gamma_lut(std::shared_ptr<const TransferFunction> fn);
-
-/** @param rgb RGB data; packed RGB 16:16:16, 48bpp, 16R, 16G, 16B,
- *  with the 2-byte value for each R/G/B component stored as
- *  little-endian; i.e. AV_PIX_FMT_RGB48LE.
- *  @param dst Buffer to fill with packed 16-bit XYZ data, i.e. first
- *  16-bit word is X, second is Y etc.
- *  @param size size of RGB image in pixels.
- *  @param size stride of RGB data in pixels.
- */
-extern void rgb_to_xyz (
-	uint8_t const * rgb,
-	uint16_t* dst,
-	dcp::Size size,
-	int stride,
-	ColourConversion const& conversion
-	);
-
-/** @param rgb RGB data; packed RGB 16:16:16, 48bpp, 16R, 16G, 16B,
- *  with the 2-byte value for each R/G/B component stored as
- *  little-endian; i.e. AV_PIX_FMT_RGB48LE.
- *  @param size size of RGB image in pixels.
- *  @param size stride of RGB data in pixels.
- */
-extern std::shared_ptr<OpenJPEGImage> rgb_to_xyz (
-	uint8_t const * rgb,
-	dcp::Size size,
-	int stride,
-	ColourConversion const& conversion
-	);
-
-
-/** @param conversion Colour conversion.
- *  @param matrix Filled in with the product of the RGB to XYZ matrix, the Bradford transform and the DCI companding.
- */
-extern void combined_rgb_to_xyz (ColourConversion const & conversion, double* matrix);
-
-}
diff --git a/src/stereo_j2k_picture_frame.cc b/src/stereo_j2k_picture_frame.cc
index 9ef91c5c..94aad06a 100644
--- a/src/stereo_j2k_picture_frame.cc
+++ b/src/stereo_j2k_picture_frame.cc
@@ -42,7 +42,6 @@
 #include "crypto_context.h"
 #include "exceptions.h"
 #include "j2k_transcode.h"
-#include "rgb_xyz.h"
 #include "stereo_j2k_picture_frame.h"
 #include "util.h"
 #include <asdcp/AS_DCP.h>
diff --git a/src/wscript b/src/wscript
index 3dc2ec12..a91e97ff 100644
--- a/src/wscript
+++ b/src/wscript
@@ -108,7 +108,6 @@ def build(bld):
              reel_stereo_picture_asset.cc
              reel_text_asset.cc
              ref.cc
-             rgb_xyz.cc
              ruby.cc
              s_gamut3_transfer_function.cc
              search.cc
@@ -227,7 +226,6 @@ def build(bld):
               reel_stereo_picture_asset.h
               reel_text_asset.h
               ref.h
-              rgb_xyz.h
               ruby.h
               s_gamut3_transfer_function.h
               scope_guard.h
diff --git a/test/colour_conversion_test.cc b/test/colour_conversion_test.cc
index 54708753..f445d2e6 100644
--- a/test/colour_conversion_test.cc
+++ b/test/colour_conversion_test.cc
@@ -31,16 +31,36 @@
     files in the program, then also delete it here.
 */
 
-#include "gamma_transfer_function.h"
+
 #include "colour_conversion.h"
+#include "gamma_transfer_function.h"
 #include "modified_gamma_transfer_function.h"
+#include "openjpeg_image.h"
+#include "piecewise_lut.h"
+#include "stream_operators.h"
+#include <boost/bind/bind.hpp>
+#include <boost/random.hpp>
+#include <boost/scoped_array.hpp>
 #include <boost/test/unit_test.hpp>
 #include <cmath>
 
+
+using std::cout;
+using std::list;
+using std::make_shared;
+using std::max;
 using std::pow;
 using std::shared_ptr;
+using std::shared_ptr;
+using std::string;
+using boost::optional;
+using boost::scoped_array;
+#if BOOST_VERSION >= 106100
+using namespace boost::placeholders;
+#endif
 using namespace dcp;
 
+
 static void
 check_gamma (shared_ptr<const TransferFunction> tf, int bit_depth, bool inverse, float gamma)
 {
@@ -154,3 +174,203 @@ BOOST_AUTO_TEST_CASE (colour_conversion_bradford_test)
 	BOOST_CHECK_CLOSE (b(2, 1), 0.0119945, 0.1);
 	BOOST_CHECK_CLOSE (b(2, 2), 0.7785377, 0.1);
 }
+
+
+static
+void
+rgb_xyz_test_case(std::function<void (uint16_t*)> write_pixel)
+{
+	srand(0);
+	dcp::Size const size(640, 480);
+
+	scoped_array<uint8_t> rgb(new uint8_t[size.width * size.height * 6]);
+	for (int y = 0; y < size.height; ++y) {
+		uint16_t* p = reinterpret_cast<uint16_t*>(rgb.get() + y * size.width * 6);
+		for (int x = 0; x < size.width; ++x) {
+			write_pixel(p);
+			p += 3;
+		}
+	}
+
+	auto xyz = dcp::rgb_to_xyz(rgb.get(), size, size.width * 6, dcp::ColourConversion::srgb_to_dcp());
+
+	for (int y = 0; y < size.height; ++y) {
+		uint16_t* p = reinterpret_cast<uint16_t*>(rgb.get() + y * size.width * 6);
+		for (int x = 0; x < size.width; ++x) {
+
+			double cr = *p++ / 65535.0;
+			double cg = *p++ / 65535.0;
+			double cb = *p++ / 65535.0;
+
+			/* Input gamma */
+
+			if (cr < 0.04045) {
+				cr /= 12.92;
+			} else {
+				cr = pow ((cr + 0.055) / 1.055, 2.4);
+			}
+
+			if (cg < 0.04045) {
+				cg /= 12.92;
+			} else {
+				cg = pow ((cg + 0.055) / 1.055, 2.4);
+			}
+
+			if (cb < 0.04045) {
+				cb /= 12.92;
+			} else {
+				cb = pow ((cb + 0.055) / 1.055, 2.4);
+			}
+
+			/* Matrix */
+
+			double cx = cr * 0.4124564 + cg * 0.3575761 + cb * 0.1804375;
+			double cy = cr * 0.2126729 + cg * 0.7151522 + cb * 0.0721750;
+			double cz = cr * 0.0193339 + cg * 0.1191920 + cb * 0.9503041;
+
+			/* Compand */
+
+			cx *= 48 / 52.37;
+			cy *= 48 / 52.37;
+			cz *= 48 / 52.37;
+
+			/* Output gamma */
+
+			cx = pow (cx, 1 / 2.6);
+			cy = pow (cy, 1 / 2.6);
+			cz = pow (cz, 1 / 2.6);
+
+			BOOST_REQUIRE_CLOSE (cx * 4095, xyz->data(0)[y * size.width + x], 1);
+			BOOST_REQUIRE_CLOSE (cy * 4095, xyz->data(1)[y * size.width + x], 1);
+			BOOST_REQUIRE_CLOSE (cz * 4095, xyz->data(2)[y * size.width + x], 1);
+		}
+	}
+}
+
+
+/** Convert a test image from sRGB to XYZ and check that the transforms are right */
+BOOST_AUTO_TEST_CASE(rgb_xyz_test)
+{
+	{
+		int counter = 0;
+		rgb_xyz_test_case([&counter](uint16_t* p) {
+			p[0] = p[1] = p[2] = (counter << 4);
+			++counter;
+			if (counter >= 4096) {
+				counter = 0;
+			}
+		});
+	}
+
+	boost::random::mt19937 rng(1);
+	boost::random::uniform_int_distribution<> dist(0, 4095);
+
+	rgb_xyz_test_case([&rng, &dist](uint16_t* p) {
+			  p[0] = dist(rng) << 4;
+			  p[1] = dist(rng) << 4;
+			  p[2] = dist(rng) << 4;
+			  });
+}
+
+
+/** Check the piecewise LUT that is used for inverse gamma calculation */
+BOOST_AUTO_TEST_CASE (rgb_xyz_lut_test)
+{
+	auto conversion = dcp::ColourConversion::rec709_to_dcp();
+	auto lut = dcp::make_inverse_gamma_lut(conversion.out_j2k());
+
+	for (double x = 0; x < 1; x += 0.000001) {
+		BOOST_CHECK(std::abs(lut.lookup(x) - lrint(pow(x, 1 / 2.6) * 4095)) < 2);
+	}
+}
+
+
+static list<string> notes;
+
+static void
+note_handler (dcp::NoteType n, string s)
+{
+	BOOST_REQUIRE_EQUAL (n, dcp::NoteType::NOTE);
+	notes.push_back (s);
+}
+
+/** Check that xyz_to_rgb clamps XYZ values correctly */
+BOOST_AUTO_TEST_CASE (xyz_rgb_range_test)
+{
+	auto xyz = make_shared<dcp::OpenJPEGImage>(dcp::Size(2, 2));
+
+	xyz->data(0)[0] = -4;
+	xyz->data(0)[1] = 6901;
+	xyz->data(0)[2] = 0;
+	xyz->data(0)[3] = 4095;
+	xyz->data(1)[0] = -4;
+	xyz->data(1)[1] = 6901;
+	xyz->data(1)[2] = 0;
+	xyz->data(1)[3] = 4095;
+	xyz->data(2)[0] = -4;
+	xyz->data(2)[1] = 6901;
+	xyz->data(2)[2] = 0;
+	xyz->data(2)[3] = 4095;
+
+	scoped_array<uint8_t> rgb (new uint8_t[2 * 2 * 6]);
+
+	notes.clear ();
+	dcp::xyz_to_rgb (
+		xyz, dcp::ColourConversion::srgb_to_dcp(), 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(), 6U);
+	auto i = notes.begin ();
+	BOOST_REQUIRE_EQUAL (*i++, "XYZ value -4 out of range");
+	BOOST_REQUIRE_EQUAL (*i++, "XYZ value -4 out of range");
+	BOOST_REQUIRE_EQUAL (*i++, "XYZ value -4 out of range");
+	BOOST_REQUIRE_EQUAL (*i++, "XYZ value 6901 out of range");
+	BOOST_REQUIRE_EQUAL (*i++, "XYZ value 6901 out of range");
+	BOOST_REQUIRE_EQUAL (*i++, "XYZ value 6901 out of range");
+
+	/* And those samples should have been clamped, so check that they give the same result
+	   as inputs at the extremes (0 and 4095).
+	*/
+
+	auto buffer = reinterpret_cast<uint16_t*> (rgb.get ());
+	BOOST_REQUIRE_EQUAL (buffer[0 * 3 + 0], buffer[2 * 3 + 1]);
+	BOOST_REQUIRE_EQUAL (buffer[0 * 3 + 1], buffer[2 * 3 + 1]);
+	BOOST_REQUIRE_EQUAL (buffer[0 * 3 + 2], buffer[2 * 3 + 2]);
+
+	BOOST_REQUIRE_EQUAL (buffer[1 * 3 + 0], buffer[3 * 3 + 0]);
+	BOOST_REQUIRE_EQUAL (buffer[1 * 3 + 1], buffer[3 * 3 + 1]);
+	BOOST_REQUIRE_EQUAL (buffer[1 * 3 + 2], buffer[3 * 3 + 2]);
+}
+
+/** Convert an image from RGB to XYZ and back again */
+BOOST_AUTO_TEST_CASE (rgb_xyz_round_trip_test)
+{
+	srand (0);
+	dcp::Size const size (640, 480);
+
+	scoped_array<uint8_t> rgb (new uint8_t[size.width * size.height * 6]);
+	for (int y = 0; y < size.height; ++y) {
+		uint16_t* p = reinterpret_cast<uint16_t*> (rgb.get() + y * size.width * 6);
+		for (int x = 0; x < size.width; ++x) {
+			/* Write a 12-bit random number for each component */
+			for (int c = 0; c < 3; ++c) {
+				*p = (rand () & 0xfff) << 4;
+				++p;
+			}
+		}
+	}
+
+	auto xyz = dcp::rgb_to_xyz(rgb.get(), size, size.width * 6, dcp::ColourConversion::srgb_to_dcp());
+	scoped_array<uint8_t> back (new uint8_t[size.width * size.height * 6]);
+	dcp::xyz_to_rgb(xyz, dcp::ColourConversion::srgb_to_dcp(), back.get(), size.width * 6);
+
+#if 0
+	uint16_t* p = reinterpret_cast<uint16_t*> (rgb.get ());
+	uint16_t* q = reinterpret_cast<uint16_t*> (back.get ());
+	for (int i = 0; i < (size.width * size.height); ++i) {
+		/* XXX: doesn't quite work */
+		// BOOST_REQUIRE_EQUAL (*p++, *q++);
+	}
+#endif
+}
diff --git a/test/decryption_test.cc b/test/decryption_test.cc
index 66402e24..061b50c8 100644
--- a/test/decryption_test.cc
+++ b/test/decryption_test.cc
@@ -47,7 +47,6 @@
 #include "reel_picture_asset.h"
 #include "reel_sound_asset.h"
 #include "reel_smpte_text_asset.h"
-#include "rgb_xyz.h"
 #include "smpte_text_asset.h"
 #include "sound_asset.h"
 #include "sound_asset_writer.h"
diff --git a/test/rgb_xyz_test.cc b/test/rgb_xyz_test.cc
deleted file mode 100644
index 59ace3d2..00000000
--- a/test/rgb_xyz_test.cc
+++ /dev/null
@@ -1,256 +0,0 @@
-/*
-    Copyright (C) 2014-2021 Carl Hetherington <cth@carlh.net>
-
-    This file is part of libdcp.
-
-    libdcp 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.
-
-    libdcp 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 libdcp.  If not, see <http://www.gnu.org/licenses/>.
-
-    In addition, as a special exception, the copyright holders give
-    permission to link the code of portions of this program with the
-    OpenSSL library under certain conditions as described in each
-    individual source file, and distribute linked combinations
-    including the two.
-
-    You must obey the GNU General Public License in all respects
-    for all of the code used other than OpenSSL.  If you modify
-    file(s) with this exception, you may extend this exception to your
-    version of the file(s), but you are not obligated to do so.  If you
-    do not wish to do so, delete this exception statement from your
-    version.  If you delete this exception statement from all source
-    files in the program, then also delete it here.
-*/
-
-
-#include "colour_conversion.h"
-#include "openjpeg_image.h"
-#include "piecewise_lut.h"
-#include "rgb_xyz.h"
-#include "stream_operators.h"
-#include <boost/bind/bind.hpp>
-#include <boost/random.hpp>
-#include <boost/scoped_array.hpp>
-#include <boost/test/unit_test.hpp>
-
-
-using std::cout;
-using std::list;
-using std::make_shared;
-using std::max;
-using std::shared_ptr;
-using std::string;
-using boost::optional;
-using boost::scoped_array;
-#if BOOST_VERSION >= 106100
-using namespace boost::placeholders;
-#endif
-
-
-static
-void
-rgb_xyz_test_case (std::function<void (uint16_t*)> write_pixel)
-{
-	srand (0);
-	dcp::Size const size (640, 480);
-
-	scoped_array<uint8_t> rgb (new uint8_t[size.width * size.height * 6]);
-	for (int y = 0; y < size.height; ++y) {
-		uint16_t* p = reinterpret_cast<uint16_t*> (rgb.get() + y * size.width * 6);
-		for (int x = 0; x < size.width; ++x) {
-			write_pixel (p);
-			p += 3;
-		}
-	}
-
-	auto xyz = dcp::rgb_to_xyz(rgb.get(), size, size.width * 6, dcp::ColourConversion::srgb_to_dcp());
-
-	for (int y = 0; y < size.height; ++y) {
-		uint16_t* p = reinterpret_cast<uint16_t*> (rgb.get() + y * size.width * 6);
-		for (int x = 0; x < size.width; ++x) {
-
-			double cr = *p++ / 65535.0;
-			double cg = *p++ / 65535.0;
-			double cb = *p++ / 65535.0;
-
-			/* Input gamma */
-
-			if (cr < 0.04045) {
-				cr /= 12.92;
-			} else {
-				cr = pow ((cr + 0.055) / 1.055, 2.4);
-			}
-
-			if (cg < 0.04045) {
-				cg /= 12.92;
-			} else {
-				cg = pow ((cg + 0.055) / 1.055, 2.4);
-			}
-
-			if (cb < 0.04045) {
-				cb /= 12.92;
-			} else {
-				cb = pow ((cb + 0.055) / 1.055, 2.4);
-			}
-
-			/* Matrix */
-
-			double cx = cr * 0.4124564 + cg * 0.3575761 + cb * 0.1804375;
-			double cy = cr * 0.2126729 + cg * 0.7151522 + cb * 0.0721750;
-			double cz = cr * 0.0193339 + cg * 0.1191920 + cb * 0.9503041;
-
-			/* Compand */
-
-			cx *= 48 / 52.37;
-			cy *= 48 / 52.37;
-			cz *= 48 / 52.37;
-
-			/* Output gamma */
-
-			cx = pow (cx, 1 / 2.6);
-			cy = pow (cy, 1 / 2.6);
-			cz = pow (cz, 1 / 2.6);
-
-			BOOST_REQUIRE_CLOSE (cx * 4095, xyz->data(0)[y * size.width + x], 1);
-			BOOST_REQUIRE_CLOSE (cy * 4095, xyz->data(1)[y * size.width + x], 1);
-			BOOST_REQUIRE_CLOSE (cz * 4095, xyz->data(2)[y * size.width + x], 1);
-		}
-	}
-}
-
-
-/** Convert a test image from sRGB to XYZ and check that the transforms are right */
-BOOST_AUTO_TEST_CASE (rgb_xyz_test)
-{
-	{
-		int counter = 0;
-		rgb_xyz_test_case ([&counter](uint16_t* p) {
-			p[0] = p[1] = p[2] = (counter << 4);
-			++counter;
-			if (counter >= 4096) {
-				counter = 0;
-			}
-		});
-	}
-
-	boost::random::mt19937 rng(1);
-	boost::random::uniform_int_distribution<> dist(0, 4095);
-
-	rgb_xyz_test_case ([&rng, &dist](uint16_t* p) {
-			   p[0] = dist(rng) << 4;
-			   p[1] = dist(rng) << 4;
-			   p[2] = dist(rng) << 4;
-			   });
-}
-
-
-/** Check the piecewise LUT that is used for inverse gamma calculation */
-BOOST_AUTO_TEST_CASE (rgb_xyz_lut_test)
-{
-	auto conversion = dcp::ColourConversion::rec709_to_dcp();
-	auto lut = dcp::make_inverse_gamma_lut(conversion.out_j2k());
-
-	for (double x = 0; x < 1; x += 0.000001) {
-		BOOST_CHECK(std::abs(lut.lookup(x) - lrint(pow(x, 1 / 2.6) * 4095)) < 2);
-	}
-}
-
-
-static list<string> notes;
-
-static void
-note_handler (dcp::NoteType n, string s)
-{
-	BOOST_REQUIRE_EQUAL (n, dcp::NoteType::NOTE);
-	notes.push_back (s);
-}
-
-/** Check that xyz_to_rgb clamps XYZ values correctly */
-BOOST_AUTO_TEST_CASE (xyz_rgb_range_test)
-{
-	auto xyz = make_shared<dcp::OpenJPEGImage>(dcp::Size(2, 2));
-
-	xyz->data(0)[0] = -4;
-	xyz->data(0)[1] = 6901;
-	xyz->data(0)[2] = 0;
-	xyz->data(0)[3] = 4095;
-	xyz->data(1)[0] = -4;
-	xyz->data(1)[1] = 6901;
-	xyz->data(1)[2] = 0;
-	xyz->data(1)[3] = 4095;
-	xyz->data(2)[0] = -4;
-	xyz->data(2)[1] = 6901;
-	xyz->data(2)[2] = 0;
-	xyz->data(2)[3] = 4095;
-
-	scoped_array<uint8_t> rgb (new uint8_t[2 * 2 * 6]);
-
-	notes.clear ();
-	dcp::xyz_to_rgb (
-		xyz, dcp::ColourConversion::srgb_to_dcp(), 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(), 6U);
-	auto i = notes.begin ();
-	BOOST_REQUIRE_EQUAL (*i++, "XYZ value -4 out of range");
-	BOOST_REQUIRE_EQUAL (*i++, "XYZ value -4 out of range");
-	BOOST_REQUIRE_EQUAL (*i++, "XYZ value -4 out of range");
-	BOOST_REQUIRE_EQUAL (*i++, "XYZ value 6901 out of range");
-	BOOST_REQUIRE_EQUAL (*i++, "XYZ value 6901 out of range");
-	BOOST_REQUIRE_EQUAL (*i++, "XYZ value 6901 out of range");
-
-	/* And those samples should have been clamped, so check that they give the same result
-	   as inputs at the extremes (0 and 4095).
-	*/
-
-	auto buffer = reinterpret_cast<uint16_t*> (rgb.get ());
-	BOOST_REQUIRE_EQUAL (buffer[0 * 3 + 0], buffer[2 * 3 + 1]);
-	BOOST_REQUIRE_EQUAL (buffer[0 * 3 + 1], buffer[2 * 3 + 1]);
-	BOOST_REQUIRE_EQUAL (buffer[0 * 3 + 2], buffer[2 * 3 + 2]);
-
-	BOOST_REQUIRE_EQUAL (buffer[1 * 3 + 0], buffer[3 * 3 + 0]);
-	BOOST_REQUIRE_EQUAL (buffer[1 * 3 + 1], buffer[3 * 3 + 1]);
-	BOOST_REQUIRE_EQUAL (buffer[1 * 3 + 2], buffer[3 * 3 + 2]);
-}
-
-/** Convert an image from RGB to XYZ and back again */
-BOOST_AUTO_TEST_CASE (rgb_xyz_round_trip_test)
-{
-	srand (0);
-	dcp::Size const size (640, 480);
-
-	scoped_array<uint8_t> rgb (new uint8_t[size.width * size.height * 6]);
-	for (int y = 0; y < size.height; ++y) {
-		uint16_t* p = reinterpret_cast<uint16_t*> (rgb.get() + y * size.width * 6);
-		for (int x = 0; x < size.width; ++x) {
-			/* Write a 12-bit random number for each component */
-			for (int c = 0; c < 3; ++c) {
-				*p = (rand () & 0xfff) << 4;
-				++p;
-			}
-		}
-	}
-
-	auto xyz = dcp::rgb_to_xyz(rgb.get(), size, size.width * 6, dcp::ColourConversion::srgb_to_dcp());
-	scoped_array<uint8_t> back (new uint8_t[size.width * size.height * 6]);
-	dcp::xyz_to_rgb(xyz, dcp::ColourConversion::srgb_to_dcp(), back.get(), size.width * 6);
-
-#if 0
-	uint16_t* p = reinterpret_cast<uint16_t*> (rgb.get ());
-	uint16_t* q = reinterpret_cast<uint16_t*> (back.get ());
-	for (int i = 0; i < (size.width * size.height); ++i) {
-		/* XXX: doesn't quite work */
-		// BOOST_REQUIRE_EQUAL (*p++, *q++);
-	}
-#endif
-}
diff --git a/test/round_trip_test.cc b/test/round_trip_test.cc
index 1e87912a..c512c6a8 100644
--- a/test/round_trip_test.cc
+++ b/test/round_trip_test.cc
@@ -31,24 +31,23 @@
     files in the program, then also delete it here.
 */
 
+
 #include "certificate.h"
+#include "certificate_chain.h"
+#include "colour_conversion.h"
+#include "cpl.h"
 #include "decrypted_kdm.h"
 #include "encrypted_kdm.h"
-#include "certificate_chain.h"
 #include "mono_j2k_picture_asset.h"
-#include "sound_asset.h"
-#include "reel.h"
-#include "test.h"
-#include "cpl.h"
-#include "mono_j2k_picture_frame.h"
-#include "certificate_chain.h"
-#include "mono_j2k_picture_asset_writer.h"
 #include "mono_j2k_picture_asset_reader.h"
-#include "reel_picture_asset.h"
-#include "reel_mono_picture_asset.h"
+#include "mono_j2k_picture_asset_writer.h"
+#include "mono_j2k_picture_frame.h"
 #include "openjpeg_image.h"
-#include "rgb_xyz.h"
-#include "colour_conversion.h"
+#include "reel.h"
+#include "reel_mono_picture_asset.h"
+#include "reel_picture_asset.h"
+#include "sound_asset.h"
+#include "test.h"
 #include <boost/test/unit_test.hpp>
 #include <boost/scoped_array.hpp>
 #include <iostream>
diff --git a/test/wscript b/test/wscript
index 5e94b1f7..ebb3a1e3 100644
--- a/test/wscript
+++ b/test/wscript
@@ -107,7 +107,6 @@ def build(bld):
                  read_change_write_test.cc
                  reel_asset_test.cc
                  recovery_test.cc
-                 rgb_xyz_test.cc
                  round_trip_test.cc
                  scope_guard_test.cc
                  search_test.cc
diff --git a/tools/dcpdumpimage.cc b/tools/dcpdumpimage.cc
index 07c32861..eeba9cc5 100644
--- a/tools/dcpdumpimage.cc
+++ b/tools/dcpdumpimage.cc
@@ -40,7 +40,6 @@
 #include "reel_picture_asset.h"
 #include "raw_convert.h"
 #include "reel.h"
-#include "rgb_xyz.h"
 #include "util.h"
 #include "warnings.h"
 LIBDCP_DISABLE_WARNINGS