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 /** @file  rgb_xyz.cc
  36  *  @brief Conversion between RGB and XYZ
  37  */
  38
  39
  40 #include "colour_conversion.h"
  41 #include "compose.hpp"
  42 #include "dcp_assert.h"
  43 #include "openjpeg_image.h"
  44 #include "rgb_xyz.h"
  45 #include "transfer_function.h"
  46 #include <cmath>
  47
  48
  49 using std::cout;
  50 using std::make_shared;
  51 using std::max;
  52 using std::min;
  53 using std::shared_ptr;
  54 using boost::optional;
  55 using namespace dcp;
  56
  57
  58 static auto constexpr DCI_COEFFICIENT = 48.0 / 52.37;
  59
  60
  61 void
  62 dcp::xyz_to_rgba (
  63         std::shared_ptr<const OpenJPEGImage> xyz_image,
  64         ColourConversion const & conversion,
  65         uint8_t* argb,
  66         int stride
  67         )
  68 {
  69         int const max_colour = pow (2, 16) - 1;
  70
  71         struct {
  72                 double x, y, z;
  73         } s;
  74
  75         struct {
  76                 double r, g, b;
  77         } d;
  78
  79         int* xyz_x = xyz_image->data (0);
  80         int* xyz_y = xyz_image->data (1);
  81         int* xyz_z = xyz_image->data (2);
  82
  83         double const * lut_in = conversion.out()->lut (12, false);
  84         double const * lut_out = conversion.in()->lut (16, true);
  85         boost::numeric::ublas::matrix<double> const matrix = conversion.xyz_to_rgb ();
  86
  87         double fast_matrix[9] = {
  88                 matrix (0, 0), matrix (0, 1), matrix (0, 2),
  89                 matrix (1, 0), matrix (1, 1), matrix (1, 2),
  90                 matrix (2, 0), matrix (2, 1), matrix (2, 2)
  91         };
  92
  93         int const height = xyz_image->size().height;
  94         int const width = xyz_image->size().width;
  95
  96         for (int y = 0; y < height; ++y) {
  97                 uint8_t* argb_line = argb;
  98                 for (int x = 0; x < width; ++x) {
  99
 100                         DCP_ASSERT (*xyz_x >= 0 && *xyz_y >= 0 && *xyz_z >= 0 && *xyz_x < 4096 && *xyz_y < 4096 && *xyz_z < 4096);
 101
 102                         /* In gamma LUT */
 103                         s.x = lut_in[*xyz_x++];
 104                         s.y = lut_in[*xyz_y++];
 105                         s.z = lut_in[*xyz_z++];
 106
 107                         /* DCI companding */
 108                         s.x /= DCI_COEFFICIENT;
 109                         s.y /= DCI_COEFFICIENT;
 110                         s.z /= DCI_COEFFICIENT;
 111
 112                         /* XYZ to RGB */
 113                         d.r = ((s.x * fast_matrix[0]) + (s.y * fast_matrix[1]) + (s.z * fast_matrix[2]));
 114                         d.g = ((s.x * fast_matrix[3]) + (s.y * fast_matrix[4]) + (s.z * fast_matrix[5]));
 115                         d.b = ((s.x * fast_matrix[6]) + (s.y * fast_matrix[7]) + (s.z * fast_matrix[8]));
 116
 117                         d.r = min (d.r, 1.0);
 118                         d.r = max (d.r, 0.0);
 119
 120                         d.g = min (d.g, 1.0);
 121                         d.g = max (d.g, 0.0);
 122
 123                         d.b = min (d.b, 1.0);
 124                         d.b = max (d.b, 0.0);
 125
 126                         /* Out gamma LUT */
 127                         *argb_line++ = lut_out[lrint(d.b * max_colour)] * 0xff;
 128                         *argb_line++ = lut_out[lrint(d.g * max_colour)] * 0xff;
 129                         *argb_line++ = lut_out[lrint(d.r * max_colour)] * 0xff;
 130                         *argb_line++ = 0xff;
 131                 }
 132
 133                 argb += stride;
 134         }
 135 }
 136
 137
 138 void
 139 dcp::xyz_to_rgb (
 140         shared_ptr<const OpenJPEGImage> xyz_image,
 141         ColourConversion const & conversion,
 142         uint8_t* rgb,
 143         int stride,
 144         optional<NoteHandler> note
 145         )
 146 {
 147         struct {
 148                 double x, y, z;
 149         } s;
 150
 151         struct {
 152                 double r, g, b;
 153         } d;
 154
 155         /* These should be 12-bit values from 0-4095 */
 156         int* xyz_x = xyz_image->data (0);
 157         int* xyz_y = xyz_image->data (1);
 158         int* xyz_z = xyz_image->data (2);
 159
 160         double const * lut_in = conversion.out()->lut (12, false);
 161         double const * lut_out = conversion.in()->lut (16, true);
 162         auto const matrix = conversion.xyz_to_rgb ();
 163
 164         double fast_matrix[9] = {
 165                 matrix (0, 0), matrix (0, 1), matrix (0, 2),
 166                 matrix (1, 0), matrix (1, 1), matrix (1, 2),
 167                 matrix (2, 0), matrix (2, 1), matrix (2, 2)
 168         };
 169
 170         int const height = xyz_image->size().height;
 171         int const width = xyz_image->size().width;
 172
 173         for (int y = 0; y < height; ++y) {
 174                 auto rgb_line = reinterpret_cast<uint16_t*> (rgb + y * stride);
 175                 for (int x = 0; x < width; ++x) {
 176
 177                         int cx = *xyz_x++;
 178                         int cy = *xyz_y++;
 179                         int cz = *xyz_z++;
 180
 181                         if (cx < 0 || cx > 4095) {
 182                                 if (note) {
 183                                         note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cx));
 184                                 }
 185                                 cx = max (min (cx, 4095), 0);
 186                         }
 187
 188                         if (cy < 0 || cy > 4095) {
 189                                 if (note) {
 190                                         note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cy));
 191                                 }
 192                                 cy = max (min (cy, 4095), 0);
 193                         }
 194
 195                         if (cz < 0 || cz > 4095) {
 196                                 if (note) {
 197                                         note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cz));
 198                                 }
 199                                 cz = max (min (cz, 4095), 0);
 200                         }
 201
 202                         /* In gamma LUT */
 203                         s.x = lut_in[cx];
 204                         s.y = lut_in[cy];
 205                         s.z = lut_in[cz];
 206
 207                         /* DCI companding */
 208                         s.x /= DCI_COEFFICIENT;
 209                         s.y /= DCI_COEFFICIENT;
 210                         s.z /= DCI_COEFFICIENT;
 211
 212                         /* XYZ to RGB */
 213                         d.r = ((s.x * fast_matrix[0]) + (s.y * fast_matrix[1]) + (s.z * fast_matrix[2]));
 214                         d.g = ((s.x * fast_matrix[3]) + (s.y * fast_matrix[4]) + (s.z * fast_matrix[5]));
 215                         d.b = ((s.x * fast_matrix[6]) + (s.y * fast_matrix[7]) + (s.z * fast_matrix[8]));
 216
 217                         d.r = min (d.r, 1.0);
 218                         d.r = max (d.r, 0.0);
 219
 220                         d.g = min (d.g, 1.0);
 221                         d.g = max (d.g, 0.0);
 222
 223                         d.b = min (d.b, 1.0);
 224                         d.b = max (d.b, 0.0);
 225
 226                         *rgb_line++ = lrint(lut_out[lrint(d.r * 65535)] * 65535);
 227                         *rgb_line++ = lrint(lut_out[lrint(d.g * 65535)] * 65535);
 228                         *rgb_line++ = lrint(lut_out[lrint(d.b * 65535)] * 65535);
 229                 }
 230         }
 231 }
 232
 233 void
 234 dcp::combined_rgb_to_xyz (ColourConversion const & conversion, double* matrix)
 235 {
 236         auto const rgb_to_xyz = conversion.rgb_to_xyz ();
 237         auto const bradford = conversion.bradford ();
 238
 239         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))
 240                 * DCI_COEFFICIENT * 65535;
 241         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))
 242                 * DCI_COEFFICIENT * 65535;
 243         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))
 244                 * DCI_COEFFICIENT * 65535;
 245         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))
 246                 * DCI_COEFFICIENT * 65535;
 247         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))
 248                 * DCI_COEFFICIENT * 65535;
 249         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))
 250                 * DCI_COEFFICIENT * 65535;
 251         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))
 252                 * DCI_COEFFICIENT * 65535;
 253         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))
 254                 * DCI_COEFFICIENT * 65535;
 255         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))
 256                 * DCI_COEFFICIENT * 65535;
 257 }
 258
 259
 260 shared_ptr<dcp::OpenJPEGImage>
 261 dcp::rgb_to_xyz (
 262         uint8_t const * rgb,
 263         dcp::Size size,
 264         int stride,
 265         ColourConversion const & conversion,
 266         optional<NoteHandler> note
 267         )
 268 {
 269         auto xyz = make_shared<OpenJPEGImage>(size);
 270
 271         struct {
 272                 double r, g, b;
 273         } s;
 274
 275         struct {
 276                 double x, y, z;
 277         } d;
 278
 279         auto const * lut_in = conversion.in()->lut (12, false);
 280         auto const * lut_out = conversion.out()->lut (16, true);
 281
 282         /* This is is the product of the RGB to XYZ matrix, the Bradford transform and the DCI companding */
 283         double fast_matrix[9];
 284         combined_rgb_to_xyz (conversion, fast_matrix);
 285
 286         int clamped = 0;
 287         int* xyz_x = xyz->data (0);
 288         int* xyz_y = xyz->data (1);
 289         int* xyz_z = xyz->data (2);
 290         for (int y = 0; y < size.height; ++y) {
 291                 auto p = reinterpret_cast<uint16_t const *> (rgb + y * stride);
 292                 for (int x = 0; x < size.width; ++x) {
 293
 294                         /* In gamma LUT (converting 16-bit to 12-bit) */
 295                         s.r = lut_in[*p++ >> 4];
 296                         s.g = lut_in[*p++ >> 4];
 297                         s.b = lut_in[*p++ >> 4];
 298
 299                         /* RGB to XYZ, Bradford transform and DCI companding */
 300                         d.x = s.r * fast_matrix[0] + s.g * fast_matrix[1] + s.b * fast_matrix[2];
 301                         d.y = s.r * fast_matrix[3] + s.g * fast_matrix[4] + s.b * fast_matrix[5];
 302                         d.z = s.r * fast_matrix[6] + s.g * fast_matrix[7] + s.b * fast_matrix[8];
 303
 304                         /* Clamp */
 305
 306                         if (d.x < 0 || d.y < 0 || d.z < 0 || d.x > 65535 || d.y > 65535 || d.z > 65535) {
 307                                 ++clamped;
 308                         }
 309
 310                         d.x = max (0.0, d.x);
 311                         d.y = max (0.0, d.y);
 312                         d.z = max (0.0, d.z);
 313                         d.x = min (65535.0, d.x);
 314                         d.y = min (65535.0, d.y);
 315                         d.z = min (65535.0, d.z);
 316
 317                         /* Out gamma LUT */
 318                         *xyz_x++ = lrint (lut_out[lrint(d.x)] * 4095);
 319                         *xyz_y++ = lrint (lut_out[lrint(d.y)] * 4095);
 320                         *xyz_z++ = lrint (lut_out[lrint(d.z)] * 4095);
 321                 }
 322         }
 323
 324         if (clamped && note) {
 325                 note.get()(NoteType::NOTE, String::compose("%1 XYZ value(s) clamped", clamped));
 326         }
 327
 328         return xyz;
 329 }