+
+ *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);
+ }
+ }
+}
+
+/** @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 of RGB image in pixels.
+ * @param stride of RGB data in pixels.
+ */
+shared_ptr<dcp::OpenJPEGImage>
+dcp::rgb_to_xyz (
+ uint8_t const * rgb,
+ dcp::Size size,
+ int stride,
+ ColourConversion const & conversion,
+ optional<NoteHandler> note
+ )
+{
+ shared_ptr<OpenJPEGImage> xyz (new OpenJPEGImage (size));
+
+ struct {
+ double r, g, b;
+ } s;
+
+ struct {
+ double x, y, z;
+ } d;
+
+ 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 ();
+
+ /* This is is the product of the RGB to XYZ matrix, the Bradford transform and the DCI companding */
+ double fast_matrix[9] = {
+ (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 * 65535,
+ (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 * 65535,
+ (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 * 65535,
+ (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 * 65535,
+ (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 * 65535,
+ (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 * 65535,
+ (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 * 65535,
+ (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 * 65535,
+ (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 * 65535
+ };
+
+ int clamped = 0;
+ int* xyz_x = xyz->data (0);
+ int* xyz_y = xyz->data (1);
+ int* xyz_z = xyz->data (2);
+ for (int y = 0; y < size.height; ++y) {
+ uint16_t const * 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 */
+
+ if (d.x < 0 || d.y < 0 || d.z < 0 || d.x > 65535 || d.y > 65535 || d.z > 65535) {
+ ++clamped;
+ }
+
+ d.x = max (0.0, d.x);
+ d.y = max (0.0, d.y);
+ d.z = max (0.0, d.z);
+ d.x = min (65535.0, d.x);
+ d.y = min (65535.0, d.y);
+ d.z = min (65535.0, d.z);
+