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         )
 280 {
 281         struct {
 282                 double r, g, b;
 283         } s;
 284
 285         struct {
 286                 double x, y, z;
 287         } d;
 288
 289         auto lut_in = conversion.in()->double_lut(0, 1, 12, false);
 290         auto lut_out = make_inverse_gamma_lut(conversion.out());
 291
 292         /* This is is the product of the RGB to XYZ matrix, the Bradford transform and the DCI companding */
 293         double fast_matrix[9];
 294         combined_rgb_to_xyz (conversion, fast_matrix);
 295
 296         for (int y = 0; y < size.height; ++y) {
 297                 auto p = reinterpret_cast<uint16_t const *> (rgb + y * stride);
 298                 for (int x = 0; x < size.width; ++x) {
 299
 300                         /* In gamma LUT (converting 16-bit to 12-bit) */
 301                         s.r = lut_in[*p++ >> 4];
 302                         s.g = lut_in[*p++ >> 4];
 303                         s.b = lut_in[*p++ >> 4];
 304
 305                         /* RGB to XYZ, Bradford transform and DCI companding */
 306                         d.x = s.r * fast_matrix[0] + s.g * fast_matrix[1] + s.b * fast_matrix[2];
 307                         d.y = s.r * fast_matrix[3] + s.g * fast_matrix[4] + s.b * fast_matrix[5];
 308                         d.z = s.r * fast_matrix[6] + s.g * fast_matrix[7] + s.b * fast_matrix[8];
 309
 310                         /* Clamp */
 311                         d.x = max (0.0, d.x);
 312                         d.y = max (0.0, d.y);
 313                         d.z = max (0.0, d.z);
 314                         d.x = min (1.0, d.x);
 315                         d.y = min (1.0, d.y);
 316                         d.z = min (1.0, d.z);
 317
 318                         /* Out gamma LUT */
 319                         *xyz_x++ = lut_out.lookup(d.x);
 320                         *xyz_y++ = lut_out.lookup(d.y);
 321                         *xyz_z++ = lut_out.lookup(d.z);
 322                 }
 323         }
 324 }
 325
 326
 327 shared_ptr<dcp::OpenJPEGImage>
 328 dcp::rgb_to_xyz (
 329         uint8_t const * rgb,
 330         dcp::Size size,
 331         int stride,
 332         ColourConversion const & conversion
 333         )
 334 {
 335         auto xyz = make_shared<OpenJPEGImage>(size);
 336
 337         int* xyz_x = xyz->data (0);
 338         int* xyz_y = xyz->data (1);
 339         int* xyz_z = xyz->data (2);
 340
 341         rgb_to_xyz_internal(rgb, xyz_x, xyz_y, xyz_z, size, stride, conversion);
 342
 343         return xyz;
 344 }
 345
 346
 347 void
 348 dcp::rgb_to_xyz (
 349         uint8_t const * rgb,
 350         uint16_t* dst,
 351         dcp::Size size,
 352         int stride,
 353         ColourConversion const & conversion
 354         )
 355 {
 356         rgb_to_xyz_internal(rgb, dst, dst, dst, size, stride, conversion);
 357 }