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()->double_lut(0, 1, 12, false);
  85         auto lut_out = conversion.in()->double_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()->double_lut(0, 1, 12, false);
 162         auto lut_out = conversion.in()->double_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, 4095);
 266 }
 267
 268
 269 template <class T>
 270 void
 271 rgb_to_xyz_internal(
 272         uint8_t const* rgb,
 273         T*& xyz_x,
 274         T*& xyz_y,
 275         T*& xyz_z,
 276         dcp::Size size,
 277         int stride,
 278         ColourConversion const& conversion,
 279         optional<NoteHandler> note
 280         )
 281 {
 282         struct {
 283                 double r, g, b;
 284         } s;
 285
 286         struct {
 287                 double x, y, z;
 288         } d;
 289
 290         auto lut_in = conversion.in()->double_lut(0, 1, 12, false);
 291         auto lut_out = make_inverse_gamma_lut(conversion.out());
 292
 293         /* This is is the product of the RGB to XYZ matrix, the Bradford transform and the DCI companding */
 294         double fast_matrix[9];
 295         combined_rgb_to_xyz (conversion, fast_matrix);
 296
 297         int clamped = 0;
 298         for (int y = 0; y < size.height; ++y) {
 299                 auto p = reinterpret_cast<uint16_t const *> (rgb + y * stride);
 300                 for (int x = 0; x < size.width; ++x) {
 301
 302                         /* In gamma LUT (converting 16-bit to 12-bit) */
 303                         s.r = lut_in[*p++ >> 4];
 304                         s.g = lut_in[*p++ >> 4];
 305                         s.b = lut_in[*p++ >> 4];
 306
 307                         /* RGB to XYZ, Bradford transform and DCI companding */
 308                         d.x = s.r * fast_matrix[0] + s.g * fast_matrix[1] + s.b * fast_matrix[2];
 309                         d.y = s.r * fast_matrix[3] + s.g * fast_matrix[4] + s.b * fast_matrix[5];
 310                         d.z = s.r * fast_matrix[6] + s.g * fast_matrix[7] + s.b * fast_matrix[8];
 311
 312                         /* Clamp */
 313
 314                         if (d.x < 0 || d.y < 0 || d.z < 0 || d.x > 1 || d.y > 1 || d.z > 1) {
 315                                 ++clamped;
 316                         }
 317
 318                         d.x = max (0.0, d.x);
 319                         d.y = max (0.0, d.y);
 320                         d.z = max (0.0, d.z);
 321                         d.x = min (1.0, d.x);
 322                         d.y = min (1.0, d.y);
 323                         d.z = min (1.0, d.z);
 324
 325                         /* Out gamma LUT */
 326                         *xyz_x++ = lut_out.lookup(d.x);
 327                         *xyz_y++ = lut_out.lookup(d.y);
 328                         *xyz_z++ = lut_out.lookup(d.z);
 329                 }
 330         }
 331
 332         if (clamped && note) {
 333                 note.get()(NoteType::NOTE, String::compose("%1 XYZ value(s) clamped", clamped));
 334         }
 335 }
 336
 337
 338 shared_ptr<dcp::OpenJPEGImage>
 339 dcp::rgb_to_xyz (
 340         uint8_t const * rgb,
 341         dcp::Size size,
 342         int stride,
 343         ColourConversion const & conversion,
 344         optional<NoteHandler> note
 345         )
 346 {
 347         auto xyz = make_shared<OpenJPEGImage>(size);
 348
 349         int* xyz_x = xyz->data (0);
 350         int* xyz_y = xyz->data (1);
 351         int* xyz_z = xyz->data (2);
 352
 353         rgb_to_xyz_internal(rgb, xyz_x, xyz_y, xyz_z, size, stride, conversion, note);
 354
 355         return xyz;
 356 }
 357
 358
 359 void
 360 dcp::rgb_to_xyz (
 361         uint8_t const * rgb,
 362         uint16_t* dst,
 363         dcp::Size size,
 364         int stride,
 365         ColourConversion const & conversion,
 366         optional<NoteHandler> note
 367         )
 368 {
 369         rgb_to_xyz_internal(rgb, dst, dst, dst, size, stride, conversion, note);
 370 }