I believe atmos_id is a unique identifier, so create it in the
[libdcp.git] / src / rgb_xyz.cc
index 2c67e62f6734f092885c721f8ef0ff2d1d9adfd4..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"
@@ -54,10 +67,11 @@ 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;
@@ -75,6 +89,12 @@ dcp::xyz_to_rgba (
        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;
 
@@ -95,9 +115,9 @@ 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);
@@ -109,14 +129,13 @@ dcp::xyz_to_rgba (
                        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;
        }
 }
 
@@ -155,9 +174,18 @@ dcp::xyz_to_rgb (
        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++;
@@ -195,9 +223,9 @@ 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);
@@ -208,18 +236,47 @@ dcp::xyz_to_rgb (
                        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 (
@@ -240,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) {
@@ -260,39 +317,28 @@ 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.0, e.x);
-                       e.y = max (0.0, e.y);
-                       e.z = max (0.0, e.z);
-                       e.x = min (65535.0, e.x);
-                       e.y = min (65535.0, e.y);
-                       e.z = min (65535.0, e.z);
+                       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);
 
                        /* 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;
-
-                       ++jn;
+                       *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);
                }
        }
 
@@ -302,28 +348,3 @@ dcp::rgb_to_xyz (
 
        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;
-}