Add OK note when PKL and CPL annotation texts match.
[libdcp.git] / src / rgb_xyz.cc
index 1110accbbe76442b0353ff0b37d77417e3b81968..c985410753af4095f0d7d1de157de179e7c42bac 100644 (file)
@@ -1,5 +1,5 @@
 /*
-    Copyright (C) 2013-2015 Carl Hetherington <cth@carlh.net>
+    Copyright (C) 2013-2021 Carl Hetherington <cth@carlh.net>
 
     This file is part of libdcp.
 
     files in the program, then also delete it here.
 */
 
-#include "rgb_xyz.h"
-#include "openjpeg_image.h"
+
+/** @file  rgb_xyz.cc
+ *  @brief Conversion between RGB and XYZ
+ */
+
+
 #include "colour_conversion.h"
-#include "transfer_function.h"
-#include "dcp_assert.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::min;
-using std::max;
+
 using std::cout;
-using boost::shared_ptr;
+using std::make_shared;
+using std::max;
+using std::min;
+using std::shared_ptr;
 using boost::optional;
 using namespace dcp;
 
-#define DCI_COEFFICIENT (48.0 / 52.37)
-
-/** 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.
- */
+
+static auto constexpr DCI_COEFFICIENT = 48.0 / 52.37;
+
+
 void
 dcp::xyz_to_rgba (
-       boost::shared_ptr<const OpenJPEGImage> xyz_image,
+       std::shared_ptr<const OpenJPEGImage> xyz_image,
        ColourConversion const & conversion,
-       uint8_t* argb
+       uint8_t* argb,
+       int stride
        )
 {
        int const max_colour = pow (2, 16) - 1;
@@ -84,8 +81,8 @@ dcp::xyz_to_rgba (
        int* xyz_y = xyz_image->data (1);
        int* xyz_z = xyz_image->data (2);
 
-       double const * lut_in = conversion.out()->lut (12, false);
-       double const * lut_out = conversion.in()->lut (16, true);
+       auto lut_in = conversion.out()->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] = {
@@ -134,20 +131,11 @@ dcp::xyz_to_rgba (
                        *argb_line++ = 0xff;
                }
 
-               /* 4 bytes per pixel */
-               argb += width * 4;
+               argb += 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).
- */
+
 void
 dcp::xyz_to_rgb (
        shared_ptr<const OpenJPEGImage> xyz_image,
@@ -170,9 +158,9 @@ dcp::xyz_to_rgb (
        int* xyz_y = xyz_image->data (1);
        int* xyz_z = xyz_image->data (2);
 
-       double const * lut_in = conversion.out()->lut (12, false);
-       double const * lut_out = conversion.in()->lut (16, true);
-       boost::numeric::ublas::matrix<double> const matrix = conversion.xyz_to_rgb ();
+       auto lut_in = conversion.out()->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),
@@ -184,7 +172,7 @@ dcp::xyz_to_rgb (
        int const width = xyz_image->size().width;
 
        for (int y = 0; y < height; ++y) {
-               uint16_t* rgb_line = reinterpret_cast<uint16_t*> (rgb + y * stride);
+               auto rgb_line = reinterpret_cast<uint16_t*> (rgb + y * stride);
                for (int x = 0; x < width; ++x) {
 
                        int cx = *xyz_x++;
@@ -193,21 +181,21 @@ dcp::xyz_to_rgb (
 
                        if (cx < 0 || cx > 4095) {
                                if (note) {
-                                       note.get() (DCP_NOTE, String::compose ("XYZ value %1 out of range", cx));
+                                       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() (DCP_NOTE, String::compose ("XYZ value %1 out of range", cy));
+                                       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() (DCP_NOTE, String::compose ("XYZ value %1 out of range", cz));
+                                       note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cz));
                                }
                                cz = max (min (cz, 4095), 0);
                        }
@@ -243,52 +231,53 @@ dcp::xyz_to_rgb (
        }
 }
 
-/** @param conversion Colour conversion.
- *  @param matrix Filled in with the product of the RGB to XYZ matrix, the Bradford transform and the DCI companding.
- */
 void
 dcp::combined_rgb_to_xyz (ColourConversion const & conversion, double* matrix)
 {
-       boost::numeric::ublas::matrix<double> const rgb_to_xyz = conversion.rgb_to_xyz ();
-       boost::numeric::ublas::matrix<double> const bradford = conversion.bradford ();
+       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 * 65535;
+               * 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 * 65535;
+               * 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 * 65535;
+               * 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 * 65535;
+               * 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 * 65535;
+               * 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 * 65535;
+               * 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 * 65535;
+               * 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 * 65535;
+               * 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 * 65535;
+               * DCI_COEFFICIENT;
 }
 
-/** @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.
- */
-shared_ptr<dcp::OpenJPEGImage>
-dcp::rgb_to_xyz (
-       uint8_t const * rgb,
+
+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,
-       optional<NoteHandler> note
+       ColourConversion const& conversion
        )
 {
-       shared_ptr<OpenJPEGImage> xyz (new OpenJPEGImage (size));
-
        struct {
                double r, g, b;
        } s;
@@ -297,19 +286,15 @@ dcp::rgb_to_xyz (
                double x, y, z;
        } d;
 
-       double const * lut_in = conversion.in()->lut (12, false);
-       double const * lut_out = conversion.out()->lut (16, true);
+       auto lut_in = conversion.in()->double_lut(0, 1, 12, false);
+       auto lut_out = make_inverse_gamma_lut(conversion.out());
 
        /* 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);
 
-       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);
+               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) */
@@ -323,53 +308,50 @@ dcp::rgb_to_xyz (
                        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);
+                       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++ = lrint (lut_out[lrint(d.x)] * 4095);
-                       *xyz_y++ = lrint (lut_out[lrint(d.y)] * 4095);
-                       *xyz_z++ = lrint (lut_out[lrint(d.z)] * 4095);
+                       *xyz_x++ = lut_out.lookup(d.x);
+                       *xyz_y++ = lut_out.lookup(d.y);
+                       *xyz_z++ = lut_out.lookup(d.z);
                }
        }
-
-       if (clamped && note) {
-               note.get() (DCP_NOTE, String::compose ("%1 XYZ value(s) clamped", clamped));
-       }
-
-       return xyz;
 }
 
 
-/** @param xyz_16 XYZ image data in packed 16:16:16, 48bpp, 16X, 16Y,
- *  16Z, with the 2-byte value for each X/Y/Z component stored as
- *  little-endian.
- */
 shared_ptr<dcp::OpenJPEGImage>
-dcp::xyz_to_xyz (uint8_t const * xyz_16, dcp::Size size, int stride)
+dcp::rgb_to_xyz (
+       uint8_t const * rgb,
+       dcp::Size size,
+       int stride,
+       ColourConversion const & conversion
+       )
 {
-       shared_ptr<OpenJPEGImage> xyz_12 (new OpenJPEGImage (size));
+       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;
+}
 
-       int jn = 0;
-       for (int y = 0; y < size.height; ++y) {
-               uint16_t const * p = reinterpret_cast<uint16_t const *> (xyz_16 + y * stride);
-               for (int x = 0; x < size.width; ++x) {
-                       /* Truncate 16-bit to 12-bit */
-                       xyz_12->data(0)[jn] = *p++ >> 4;
-                       xyz_12->data(1)[jn] = *p++ >> 4;
-                       xyz_12->data(2)[jn] = *p++ >> 4;
-                       ++jn;
-               }
-       }
 
-       return xyz_12;
+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);
 }