src/rgb_xyz.cc

   1 /*
   2     Copyright (C) 2013-2021 Carl Hetherington <cth@carlh.net>
   3
   4     This file is part of libdcp.
   5
   6     libdcp is free software; you can redistribute it and/or modify
   7     it under the terms of the GNU General Public License as published by
   8     the Free Software Foundation; either version 2 of the License, or
   9     (at your option) any later version.
  10
  11     libdcp is distributed in the hope that it will be useful,
  12     but WITHOUT ANY WARRANTY; without even the implied warranty of
  13     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14     GNU General Public License for more details.
  15
  16     You should have received a copy of the GNU General Public License
  17     along with libdcp.  If not, see <http://www.gnu.org/licenses/>.
  18
  19     In addition, as a special exception, the copyright holders give
  20     permission to link the code of portions of this program with the
  21     OpenSSL library under certain conditions as described in each
  22     individual source file, and distribute linked combinations
  23     including the two.
  24
  25     You must obey the GNU General Public License in all respects
  26     for all of the code used other than OpenSSL.  If you modify
  27     file(s) with this exception, you may extend this exception to your
  28     version of the file(s), but you are not obligated to do so.  If you
  29     do not wish to do so, delete this exception statement from your
  30     version.  If you delete this exception statement from all source
  31     files in the program, then also delete it here.
  32 */
  33
  34
  35 #include "colour_conversion.h"
  36 #include "compose.hpp"
  37 #include "dcp_assert.h"
  38 #include "openjpeg_image.h"
  39 #include "rgb_xyz.h"
  40 #include "transfer_function.h"
  41 #include <cmath>
  42
  43
  44 using std::cout;
  45 using std::make_shared;
  46 using std::max;
  47 using std::min;
  48 using std::shared_ptr;
  49 using boost::optional;
  50 using namespace dcp;
  51
  52
  53 static auto constexpr DCI_COEFFICIENT = 48.0 / 52.37;
  54
  55
  56 void
  57 dcp::xyz_to_rgba (
  58         std::shared_ptr<const OpenJPEGImage> xyz_image,
  59         ColourConversion const & conversion,
  60         uint8_t* argb,
  61         int stride
  62         )
  63 {
  64         int const max_colour = pow (2, 16) - 1;
  65
  66         struct {
  67                 double x, y, z;
  68         } s;
  69
  70         struct {
  71                 double r, g, b;
  72         } d;
  73
  74         int* xyz_x = xyz_image->data (0);
  75         int* xyz_y = xyz_image->data (1);
  76         int* xyz_z = xyz_image->data (2);
  77
  78         double const * lut_in = conversion.out()->lut (12, false);
  79         double const * lut_out = conversion.in()->lut (16, true);
  80         boost::numeric::ublas::matrix<double> const matrix = conversion.xyz_to_rgb ();
  81
  82         double fast_matrix[9] = {
  83                 matrix (0, 0), matrix (0, 1), matrix (0, 2),
  84                 matrix (1, 0), matrix (1, 1), matrix (1, 2),
  85                 matrix (2, 0), matrix (2, 1), matrix (2, 2)
  86         };
  87
  88         int const height = xyz_image->size().height;
  89         int const width = xyz_image->size().width;
  90
  91         for (int y = 0; y < height; ++y) {
  92                 uint8_t* argb_line = argb;
  93                 for (int x = 0; x < width; ++x) {
  94
  95                         DCP_ASSERT (*xyz_x >= 0 && *xyz_y >= 0 && *xyz_z >= 0 && *xyz_x < 4096 && *xyz_y < 4096 && *xyz_z < 4096);
  96
  97                         /* In gamma LUT */
  98                         s.x = lut_in[*xyz_x++];
  99                         s.y = lut_in[*xyz_y++];
 100                         s.z = lut_in[*xyz_z++];
 101
 102                         /* DCI companding */
 103                         s.x /= DCI_COEFFICIENT;
 104                         s.y /= DCI_COEFFICIENT;
 105                         s.z /= DCI_COEFFICIENT;
 106
 107                         /* XYZ to RGB */
 108                         d.r = ((s.x * fast_matrix[0]) + (s.y * fast_matrix[1]) + (s.z * fast_matrix[2]));
 109                         d.g = ((s.x * fast_matrix[3]) + (s.y * fast_matrix[4]) + (s.z * fast_matrix[5]));
 110                         d.b = ((s.x * fast_matrix[6]) + (s.y * fast_matrix[7]) + (s.z * fast_matrix[8]));
 111
 112                         d.r = min (d.r, 1.0);
 113                         d.r = max (d.r, 0.0);
 114
 115                         d.g = min (d.g, 1.0);
 116                         d.g = max (d.g, 0.0);
 117
 118                         d.b = min (d.b, 1.0);
 119                         d.b = max (d.b, 0.0);
 120
 121                         /* Out gamma LUT */
 122                         *argb_line++ = lut_out[lrint(d.b * max_colour)] * 0xff;
 123                         *argb_line++ = lut_out[lrint(d.g * max_colour)] * 0xff;
 124                         *argb_line++ = lut_out[lrint(d.r * max_colour)] * 0xff;
 125                         *argb_line++ = 0xff;
 126                 }
 127
 128                 argb += stride;
 129         }
 130 }
 131
 132
 133 void
 134 dcp::xyz_to_rgb (
 135         shared_ptr<const OpenJPEGImage> xyz_image,
 136         ColourConversion const & conversion,
 137         uint8_t* rgb,
 138         int stride,
 139         optional<NoteHandler> note
 140         )
 141 {
 142         struct {
 143                 double x, y, z;
 144         } s;
 145
 146         struct {
 147                 double r, g, b;
 148         } d;
 149
 150         /* These should be 12-bit values from 0-4095 */
 151         int* xyz_x = xyz_image->data (0);
 152         int* xyz_y = xyz_image->data (1);
 153         int* xyz_z = xyz_image->data (2);
 154
 155         double const * lut_in = conversion.out()->lut (12, false);
 156         double const * lut_out = conversion.in()->lut (16, true);
 157         boost::numeric::ublas::matrix<double> const matrix = conversion.xyz_to_rgb ();
 158
 159         double fast_matrix[9] = {
 160                 matrix (0, 0), matrix (0, 1), matrix (0, 2),
 161                 matrix (1, 0), matrix (1, 1), matrix (1, 2),
 162                 matrix (2, 0), matrix (2, 1), matrix (2, 2)
 163         };
 164
 165         int const height = xyz_image->size().height;
 166         int const width = xyz_image->size().width;
 167
 168         for (int y = 0; y < height; ++y) {
 169                 uint16_t* rgb_line = reinterpret_cast<uint16_t*> (rgb + y * stride);
 170                 for (int x = 0; x < width; ++x) {
 171
 172                         int cx = *xyz_x++;
 173                         int cy = *xyz_y++;
 174                         int cz = *xyz_z++;
 175
 176                         if (cx < 0 || cx > 4095) {
 177                                 if (note) {
 178                                         note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cx));
 179                                 }
 180                                 cx = max (min (cx, 4095), 0);
 181                         }
 182
 183                         if (cy < 0 || cy > 4095) {
 184                                 if (note) {
 185                                         note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cy));
 186                                 }
 187                                 cy = max (min (cy, 4095), 0);
 188                         }
 189
 190                         if (cz < 0 || cz > 4095) {
 191                                 if (note) {
 192                                         note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cz));
 193                                 }
 194                                 cz = max (min (cz, 4095), 0);
 195                         }
 196
 197                         /* In gamma LUT */
 198                         s.x = lut_in[cx];
 199                         s.y = lut_in[cy];
 200                         s.z = lut_in[cz];
 201
 202                         /* DCI companding */
 203                         s.x /= DCI_COEFFICIENT;
 204                         s.y /= DCI_COEFFICIENT;
 205                         s.z /= DCI_COEFFICIENT;
 206
 207                         /* XYZ to RGB */
 208                         d.r = ((s.x * fast_matrix[0]) + (s.y * fast_matrix[1]) + (s.z * fast_matrix[2]));
 209                         d.g = ((s.x * fast_matrix[3]) + (s.y * fast_matrix[4]) + (s.z * fast_matrix[5]));
 210                         d.b = ((s.x * fast_matrix[6]) + (s.y * fast_matrix[7]) + (s.z * fast_matrix[8]));
 211
 212                         d.r = min (d.r, 1.0);
 213                         d.r = max (d.r, 0.0);
 214
 215                         d.g = min (d.g, 1.0);
 216                         d.g = max (d.g, 0.0);
 217
 218                         d.b = min (d.b, 1.0);
 219                         d.b = max (d.b, 0.0);
 220
 221                         *rgb_line++ = lrint(lut_out[lrint(d.r * 65535)] * 65535);
 222                         *rgb_line++ = lrint(lut_out[lrint(d.g * 65535)] * 65535);
 223                         *rgb_line++ = lrint(lut_out[lrint(d.b * 65535)] * 65535);
 224                 }
 225         }
 226 }
 227
 228 void
 229 dcp::combined_rgb_to_xyz (ColourConversion const & conversion, double* matrix)
 230 {
 231         boost::numeric::ublas::matrix<double> const rgb_to_xyz = conversion.rgb_to_xyz ();
 232         boost::numeric::ublas::matrix<double> const bradford = conversion.bradford ();
 233
 234         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))
 235                 * DCI_COEFFICIENT * 65535;
 236         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))
 237                 * DCI_COEFFICIENT * 65535;
 238         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))
 239                 * DCI_COEFFICIENT * 65535;
 240         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))
 241                 * DCI_COEFFICIENT * 65535;
 242         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))
 243                 * DCI_COEFFICIENT * 65535;
 244         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))
 245                 * DCI_COEFFICIENT * 65535;
 246         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))
 247                 * DCI_COEFFICIENT * 65535;
 248         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))
 249                 * DCI_COEFFICIENT * 65535;
 250         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))
 251                 * DCI_COEFFICIENT * 65535;
 252 }
 253
 254
 255 shared_ptr<dcp::OpenJPEGImage>
 256 dcp::rgb_to_xyz (
 257         uint8_t const * rgb,
 258         dcp::Size size,
 259         int stride,
 260         ColourConversion const & conversion,
 261         optional<NoteHandler> note
 262         )
 263 {
 264         auto xyz = make_shared<OpenJPEGImage>(size);
 265
 266         struct {
 267                 double r, g, b;
 268         } s;
 269
 270         struct {
 271                 double x, y, z;
 272         } d;
 273
 274         auto const * lut_in = conversion.in()->lut (12, false);
 275         auto const * lut_out = conversion.out()->lut (16, true);
 276
 277         /* This is is the product of the RGB to XYZ matrix, the Bradford transform and the DCI companding */
 278         double fast_matrix[9];
 279         combined_rgb_to_xyz (conversion, fast_matrix);
 280
 281         int clamped = 0;
 282         int* xyz_x = xyz->data (0);
 283         int* xyz_y = xyz->data (1);
 284         int* xyz_z = xyz->data (2);
 285         for (int y = 0; y < size.height; ++y) {
 286                 auto p = reinterpret_cast<uint16_t const *> (rgb + y * stride);
 287                 for (int x = 0; x < size.width; ++x) {
 288
 289                         /* In gamma LUT (converting 16-bit to 12-bit) */
 290                         s.r = lut_in[*p++ >> 4];
 291                         s.g = lut_in[*p++ >> 4];
 292                         s.b = lut_in[*p++ >> 4];
 293
 294                         /* RGB to XYZ, Bradford transform and DCI companding */
 295                         d.x = s.r * fast_matrix[0] + s.g * fast_matrix[1] + s.b * fast_matrix[2];
 296                         d.y = s.r * fast_matrix[3] + s.g * fast_matrix[4] + s.b * fast_matrix[5];
 297                         d.z = s.r * fast_matrix[6] + s.g * fast_matrix[7] + s.b * fast_matrix[8];
 298
 299                         /* Clamp */
 300
 301                         if (d.x < 0 || d.y < 0 || d.z < 0 || d.x > 65535 || d.y > 65535 || d.z > 65535) {
 302                                 ++clamped;
 303                         }
 304
 305                         d.x = max (0.0, d.x);
 306                         d.y = max (0.0, d.y);
 307                         d.z = max (0.0, d.z);
 308                         d.x = min (65535.0, d.x);
 309                         d.y = min (65535.0, d.y);
 310                         d.z = min (65535.0, d.z);
 311
 312                         /* Out gamma LUT */
 313                         *xyz_x++ = lrint (lut_out[lrint(d.x)] * 4095);
 314                         *xyz_y++ = lrint (lut_out[lrint(d.y)] * 4095);
 315                         *xyz_z++ = lrint (lut_out[lrint(d.z)] * 4095);
 316                 }
 317         }
 318
 319         if (clamped && note) {
 320                 note.get()(NoteType::NOTE, String::compose("%1 XYZ value(s) clamped", clamped));
 321         }
 322
 323         return xyz;
 324 }