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