Add LocalTime::add_minutes.
[libdcp.git] / src / rgb_xyz.cc
index 08a4e838a78867338dd0bec69b426d9cc284c0d6..d0774e7c1bcd142b002a9850a8058643b1190728 100644 (file)
@@ -1,25 +1,38 @@
 /*
-    Copyright (C) 2013-2014 Carl Hetherington <cth@carlh.net>
+    Copyright (C) 2013-2015 Carl Hetherington <cth@carlh.net>
 
-    This program is free software; you can redistribute it and/or modify
+    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.
 
-    This program is distributed in the hope that it will be useful,
+    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 this program; if not, write to the Free Software
-    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
-
+    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 "rgb_xyz.h"
 #include "openjpeg_image.h"
-#include "colour_matrix.h"
 #include "colour_conversion.h"
 #include "transfer_function.h"
 #include "dcp_assert.h"
@@ -28,6 +41,7 @@
 
 using std::min;
 using std::max;
+using std::cout;
 using boost::shared_ptr;
 using boost::optional;
 using namespace dcp;
@@ -53,36 +67,43 @@ void
 dcp::xyz_to_rgba (
        boost::shared_ptr<const OpenJPEGImage> xyz_image,
        ColourConversion const & conversion,
-       uint8_t* argb
+       uint8_t* argb,
+       int stride
        )
 {
-       int const max_colour = pow (2, 12) - 1;
+       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);
 
-       double const * lut_in = conversion.in()->lut (16, true);
-       double const * lut_out = conversion.out()->lut (12, false);
+       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 ();
 
+       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++];
@@ -94,28 +115,27 @@ dcp::xyz_to_rgba (
                        s.z /= DCI_COEFFICIENT;
 
                        /* XYZ to RGB */
-                       d.r = ((s.x * matrix(0, 0)) + (s.y * matrix(0, 1)) + (s.z * matrix(0, 2)));
-                       d.g = ((s.x * matrix(1, 0)) + (s.y * matrix(1, 1)) + (s.z * matrix(1, 2)));
-                       d.b = ((s.x * matrix(2, 0)) + (s.y * matrix(2, 1)) + (s.z * matrix(2, 2)));
-                       
+                       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[int(rint(d.b * max_colour))] * 0xff;
-                       *argb_line++ = lut_out[int(rint(d.g * max_colour))] * 0xff;
-                       *argb_line++ = lut_out[int(rint(d.r * max_colour))] * 0xff;
+                       *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;
                }
 
-               /* 4 bytes per pixel */
-               argb += width * 4;
+               argb += stride;
        }
 }
 
@@ -140,7 +160,7 @@ dcp::xyz_to_rgb (
        struct {
                double x, y, z;
        } s;
-       
+
        struct {
                double r, g, b;
        } d;
@@ -150,13 +170,22 @@ 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, true);
-       double const * lut_out = conversion.out()->lut (16, false);
+       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 ();
 
-       for (int y = 0; y < xyz_image->size().height; ++y) {
+       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) {
                uint16_t* rgb_line = reinterpret_cast<uint16_t*> (rgb + y * stride);
-               for (int x = 0; x < xyz_image->size().width; ++x) {
+               for (int x = 0; x < width; ++x) {
 
                        int cx = *xyz_x++;
                        int cy = *xyz_y++;
@@ -182,7 +211,7 @@ dcp::xyz_to_rgb (
                                }
                                cz = max (min (cz, 4095), 0);
                        }
-                       
+
                        /* In gamma LUT */
                        s.x = lut_in[cx];
                        s.y = lut_in[cy];
@@ -194,31 +223,60 @@ dcp::xyz_to_rgb (
                        s.z /= DCI_COEFFICIENT;
 
                        /* XYZ to RGB */
-                       d.r = ((s.x * matrix(0, 0)) + (s.y * matrix(0, 1)) + (s.z * matrix(0, 2)));
-                       d.g = ((s.x * matrix(1, 0)) + (s.y * matrix(1, 1)) + (s.z * matrix(1, 2)));
-                       d.b = ((s.x * matrix(2, 0)) + (s.y * matrix(2, 1)) + (s.z * matrix(2, 2)));
-                       
+                       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++ = rint(lut_out[int(rint(d.r * 65535))] * 65535);
-                       *rgb_line++ = rint(lut_out[int(rint(d.g * 65535))] * 65535);
-                       *rgb_line++ = rint(lut_out[int(rint(d.b * 65535))] * 65535);
+                       *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 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 ();
+
+       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;
+       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;
+       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;
+       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;
+       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;
+       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;
+       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;
+       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;
+       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;
+}
+
 /** @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.
+ *  @param size size of RGB image in pixels.
+ *  @param size stride of RGB data in pixels.
  */
 shared_ptr<dcp::OpenJPEGImage>
 dcp::rgb_to_xyz (
@@ -239,17 +297,17 @@ dcp::rgb_to_xyz (
                double x, y, z;
        } d;
 
-       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 ();
+
+       /* 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 jn = 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) {
@@ -259,70 +317,34 @@ dcp::rgb_to_xyz (
                        s.g = lut_in[*p++ >> 4];
                        s.b = lut_in[*p++ >> 4];
 
-                       /* RGB to XYZ Matrix */
-                       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 */
-                       e.x = e.x * DCI_COEFFICIENT * 65535;
-                       e.y = e.y * DCI_COEFFICIENT * 65535;
-                       e.z = e.z * DCI_COEFFICIENT * 65535;
+                       /* 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 (e.x < 0 || e.y < 0 || e.z < 0 || e.x > 65535 || e.y > 65535 || e.z > 65535) {
+                       if (d.x < 0 || d.y < 0 || d.z < 0 || d.x > 65535 || d.y > 65535 || d.z > 65535) {
                                ++clamped;
                        }
-                       
-                       e.x = max (0.0d, e.x);
-                       e.y = max (0.0d, e.y);
-                       e.z = max (0.0d, e.z);
-                       e.x = min (65535.0d, e.x);
-                       e.y = min (65535.0d, e.y);
-                       e.z = min (65535.0d, e.z);
 
-                       /* Out gamma LUT */
-                       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;
+                       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);
 
-                       ++jn;
+                       /* 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);
                }
        }
 
-       if (clamped) {
+       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)
-{
-       shared_ptr<OpenJPEGImage> xyz_12 (new OpenJPEGImage (size));
-
-       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;
-}