src/rgb_xyz.cc

   1 /*
   2     Copyright (C) 2013-2015 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 #include "rgb_xyz.h"
  35 #include "openjpeg_image.h"
  36 #include "colour_matrix.h"
  37 #include "colour_conversion.h"
  38 #include "transfer_function.h"
  39 #include "dcp_assert.h"
  40 #include "compose.hpp"
  41 #include <cmath>
  42
  43 using std::min;
  44 using std::max;
  45 using std::cout;
  46 using boost::shared_ptr;
  47 using boost::optional;
  48 using namespace dcp;
  49
  50 #define DCI_COEFFICIENT (48.0 / 52.37)
  51
  52 /** Convert an XYZ image to RGBA.
  53  *  @param xyz_image Image in XYZ.
  54  *  @param conversion Colour conversion to use.
  55  *  @param argb Buffer to fill with RGBA data.  The format of the data is:
  56  *
  57  *  <pre>
  58  *  Byte   /- 0 -------|- 1 --------|- 2 --------|- 3 --------|- 4 --------|- 5 --------| ...
  59  *         |(0, 0) Blue|(0, 0)Green |(0, 0) Red  |(0, 0) Alpha|(0, 1) Blue |(0, 1) Green| ...
  60  *  </pre>
  61  *
  62  *  So that the first byte is the blue component of the pixel at x=0, y=0, the second
  63  *  is the green component, and so on.
  64  *
  65  *  Lines are packed so that the second row directly follows the first.
  66  */
  67 void
  68 dcp::xyz_to_rgba (
  69         boost::shared_ptr<const OpenJPEGImage> xyz_image,
  70         ColourConversion const & conversion,
  71         uint8_t* argb
  72         )
  73 {
  74         int const max_colour = pow (2, 16) - 1;
  75
  76         struct {
  77                 double x, y, z;
  78         } s;
  79
  80         struct {
  81                 double r, g, b;
  82         } d;
  83
  84         int* xyz_x = xyz_image->data (0);
  85         int* xyz_y = xyz_image->data (1);
  86         int* xyz_z = xyz_image->data (2);
  87
  88         double const * lut_in = conversion.out()->lut (12, false);
  89         double const * lut_out = conversion.in()->lut (16, true);
  90         boost::numeric::ublas::matrix<double> const matrix = conversion.xyz_to_rgb ();
  91
  92         double fast_matrix[9] = {
  93                 matrix (0, 0), matrix (0, 1), matrix (0, 2),
  94                 matrix (1, 0), matrix (1, 1), matrix (1, 2),
  95                 matrix (2, 0), matrix (2, 1), matrix (2, 2)
  96         };
  97
  98         int const height = xyz_image->size().height;
  99         int const width = xyz_image->size().width;
 100
 101         for (int y = 0; y < height; ++y) {
 102                 uint8_t* argb_line = argb;
 103                 for (int x = 0; x < width; ++x) {
 104
 105                         DCP_ASSERT (*xyz_x >= 0 && *xyz_y >= 0 && *xyz_z >= 0 && *xyz_x < 4096 && *xyz_y < 4096 && *xyz_z < 4096);
 106
 107                         /* In gamma LUT */
 108                         s.x = lut_in[*xyz_x++];
 109                         s.y = lut_in[*xyz_y++];
 110                         s.z = lut_in[*xyz_z++];
 111
 112                         /* DCI companding */
 113                         s.x /= DCI_COEFFICIENT;
 114                         s.y /= DCI_COEFFICIENT;
 115                         s.z /= DCI_COEFFICIENT;
 116
 117                         /* XYZ to RGB */
 118                         d.r = ((s.x * fast_matrix[0]) + (s.y * fast_matrix[1]) + (s.z * fast_matrix[2]));
 119                         d.g = ((s.x * fast_matrix[3]) + (s.y * fast_matrix[4]) + (s.z * fast_matrix[5]));
 120                         d.b = ((s.x * fast_matrix[6]) + (s.y * fast_matrix[7]) + (s.z * fast_matrix[8]));
 121
 122                         d.r = min (d.r, 1.0);
 123                         d.r = max (d.r, 0.0);
 124
 125                         d.g = min (d.g, 1.0);
 126                         d.g = max (d.g, 0.0);
 127
 128                         d.b = min (d.b, 1.0);
 129                         d.b = max (d.b, 0.0);
 130
 131                         /* Out gamma LUT */
 132                         *argb_line++ = lut_out[lrint(d.b * max_colour)] * 0xff;
 133                         *argb_line++ = lut_out[lrint(d.g * max_colour)] * 0xff;
 134                         *argb_line++ = lut_out[lrint(d.r * max_colour)] * 0xff;
 135                         *argb_line++ = 0xff;
 136                 }
 137
 138                 /* 4 bytes per pixel */
 139                 argb += width * 4;
 140         }
 141 }
 142
 143 /** Convert an XYZ image to 48bpp RGB.
 144  *  @param xyz_image Frame in XYZ.
 145  *  @param conversion Colour conversion to use.
 146  *  @param rgb Buffer to fill with RGB data.  Format is packed RGB
 147  *  16:16:16, 48bpp, 16R, 16G, 16B, with the 2-byte value for each
 148  *  R/G/B component stored as little-endian; i.e. AV_PIX_FMT_RGB48LE.
 149  *  @param stride Stride for RGB data in bytes.
 150  *  @param note Optional handler for any notes that may be made during the conversion (e.g. when clamping occurs).
 151  */
 152 void
 153 dcp::xyz_to_rgb (
 154         shared_ptr<const OpenJPEGImage> xyz_image,
 155         ColourConversion const & conversion,
 156         uint8_t* rgb,
 157         int stride,
 158         optional<NoteHandler> note
 159         )
 160 {
 161         struct {
 162                 double x, y, z;
 163         } s;
 164
 165         struct {
 166                 double r, g, b;
 167         } d;
 168
 169         /* These should be 12-bit values from 0-4095 */
 170         int* xyz_x = xyz_image->data (0);
 171         int* xyz_y = xyz_image->data (1);
 172         int* xyz_z = xyz_image->data (2);
 173
 174         double const * lut_in = conversion.out()->lut (12, false);
 175         double const * lut_out = conversion.in()->lut (16, true);
 176         boost::numeric::ublas::matrix<double> const matrix = conversion.xyz_to_rgb ();
 177
 178         double fast_matrix[9] = {
 179                 matrix (0, 0), matrix (0, 1), matrix (0, 2),
 180                 matrix (1, 0), matrix (1, 1), matrix (1, 2),
 181                 matrix (2, 0), matrix (2, 1), matrix (2, 2)
 182         };
 183
 184         int const height = xyz_image->size().height;
 185         int const width = xyz_image->size().width;
 186
 187         for (int y = 0; y < height; ++y) {
 188                 uint16_t* rgb_line = reinterpret_cast<uint16_t*> (rgb + y * stride);
 189                 for (int x = 0; x < width; ++x) {
 190
 191                         int cx = *xyz_x++;
 192                         int cy = *xyz_y++;
 193                         int cz = *xyz_z++;
 194
 195                         if (cx < 0 || cx > 4095) {
 196                                 if (note) {
 197                                         note.get() (DCP_NOTE, String::compose ("XYZ value %1 out of range", cx));
 198                                 }
 199                                 cx = max (min (cx, 4095), 0);
 200                         }
 201
 202                         if (cy < 0 || cy > 4095) {
 203                                 if (note) {
 204                                         note.get() (DCP_NOTE, String::compose ("XYZ value %1 out of range", cy));
 205                                 }
 206                                 cy = max (min (cy, 4095), 0);
 207                         }
 208
 209                         if (cz < 0 || cz > 4095) {
 210                                 if (note) {
 211                                         note.get() (DCP_NOTE, String::compose ("XYZ value %1 out of range", cz));
 212                                 }
 213                                 cz = max (min (cz, 4095), 0);
 214                         }
 215
 216                         /* In gamma LUT */
 217                         s.x = lut_in[cx];
 218                         s.y = lut_in[cy];
 219                         s.z = lut_in[cz];
 220
 221                         /* DCI companding */
 222                         s.x /= DCI_COEFFICIENT;
 223                         s.y /= DCI_COEFFICIENT;
 224                         s.z /= DCI_COEFFICIENT;
 225
 226                         /* XYZ to RGB */
 227                         d.r = ((s.x * fast_matrix[0]) + (s.y * fast_matrix[1]) + (s.z * fast_matrix[2]));
 228                         d.g = ((s.x * fast_matrix[3]) + (s.y * fast_matrix[4]) + (s.z * fast_matrix[5]));
 229                         d.b = ((s.x * fast_matrix[6]) + (s.y * fast_matrix[7]) + (s.z * fast_matrix[8]));
 230
 231                         d.r = min (d.r, 1.0);
 232                         d.r = max (d.r, 0.0);
 233
 234                         d.g = min (d.g, 1.0);
 235                         d.g = max (d.g, 0.0);
 236
 237                         d.b = min (d.b, 1.0);
 238                         d.b = max (d.b, 0.0);
 239
 240                         *rgb_line++ = lrint(lut_out[lrint(d.r * 65535)] * 65535);
 241                         *rgb_line++ = lrint(lut_out[lrint(d.g * 65535)] * 65535);
 242                         *rgb_line++ = lrint(lut_out[lrint(d.b * 65535)] * 65535);
 243                 }
 244         }
 245 }
 246
 247 /** @param conversion Colour conversion.
 248  *  @param matrix Filled in with the product of the RGB to XYZ matrix, the Bradford transform and the DCI companding.
 249  */
 250 void
 251 dcp::combined_rgb_to_xyz (ColourConversion const & conversion, double* matrix)
 252 {
 253         boost::numeric::ublas::matrix<double> const rgb_to_xyz = conversion.rgb_to_xyz ();
 254         boost::numeric::ublas::matrix<double> const bradford = conversion.bradford ();
 255
 256         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))
 257                 * DCI_COEFFICIENT * 65535;
 258         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))
 259                 * DCI_COEFFICIENT * 65535;
 260         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))
 261                 * DCI_COEFFICIENT * 65535;
 262         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))
 263                 * DCI_COEFFICIENT * 65535;
 264         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))
 265                 * DCI_COEFFICIENT * 65535;
 266         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))
 267                 * DCI_COEFFICIENT * 65535;
 268         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))
 269                 * DCI_COEFFICIENT * 65535;
 270         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))
 271                 * DCI_COEFFICIENT * 65535;
 272         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))
 273                 * DCI_COEFFICIENT * 65535;
 274 }
 275
 276 /** @param rgb RGB data; packed RGB 16:16:16, 48bpp, 16R, 16G, 16B,
 277  *  with the 2-byte value for each R/G/B component stored as
 278  *  little-endian; i.e. AV_PIX_FMT_RGB48LE.
 279  *  @param size size of RGB image in pixels.
 280  *  @param size stride of RGB data in pixels.
 281  */
 282 shared_ptr<dcp::OpenJPEGImage>
 283 dcp::rgb_to_xyz (
 284         uint8_t const * rgb,
 285         dcp::Size size,
 286         int stride,
 287         ColourConversion const & conversion,
 288         optional<NoteHandler> note
 289         )
 290 {
 291         shared_ptr<OpenJPEGImage> xyz (new OpenJPEGImage (size));
 292
 293         struct {
 294                 double r, g, b;
 295         } s;
 296
 297         struct {
 298                 double x, y, z;
 299         } d;
 300
 301         double const * lut_in = conversion.in()->lut (12, false);
 302         double const * lut_out = conversion.out()->lut (16, true);
 303
 304         /* This is is the product of the RGB to XYZ matrix, the Bradford transform and the DCI companding */
 305         double fast_matrix[9];
 306         combined_rgb_to_xyz (conversion, fast_matrix);
 307
 308         int clamped = 0;
 309         int* xyz_x = xyz->data (0);
 310         int* xyz_y = xyz->data (1);
 311         int* xyz_z = xyz->data (2);
 312         for (int y = 0; y < size.height; ++y) {
 313                 uint16_t const * p = reinterpret_cast<uint16_t const *> (rgb + y * stride);
 314                 for (int x = 0; x < size.width; ++x) {
 315
 316                         /* In gamma LUT (converting 16-bit to 12-bit) */
 317                         s.r = lut_in[*p++ >> 4];
 318                         s.g = lut_in[*p++ >> 4];
 319                         s.b = lut_in[*p++ >> 4];
 320
 321                         /* RGB to XYZ, Bradford transform and DCI companding */
 322                         d.x = s.r * fast_matrix[0] + s.g * fast_matrix[1] + s.b * fast_matrix[2];
 323                         d.y = s.r * fast_matrix[3] + s.g * fast_matrix[4] + s.b * fast_matrix[5];
 324                         d.z = s.r * fast_matrix[6] + s.g * fast_matrix[7] + s.b * fast_matrix[8];
 325
 326                         /* Clamp */
 327
 328                         if (d.x < 0 || d.y < 0 || d.z < 0 || d.x > 65535 || d.y > 65535 || d.z > 65535) {
 329                                 ++clamped;
 330                         }
 331
 332                         d.x = max (0.0, d.x);
 333                         d.y = max (0.0, d.y);
 334                         d.z = max (0.0, d.z);
 335                         d.x = min (65535.0, d.x);
 336                         d.y = min (65535.0, d.y);
 337                         d.z = min (65535.0, d.z);
 338
 339                         /* Out gamma LUT */
 340                         *xyz_x++ = lrint (lut_out[lrint(d.x)] * 4095);
 341                         *xyz_y++ = lrint (lut_out[lrint(d.y)] * 4095);
 342                         *xyz_z++ = lrint (lut_out[lrint(d.z)] * 4095);
 343                 }
 344         }
 345
 346         if (clamped && note) {
 347                 note.get() (DCP_NOTE, String::compose ("%1 XYZ value(s) clamped", clamped));
 348         }
 349
 350         return xyz;
 351 }
 352
 353
 354 /** @param xyz_16 XYZ image data in packed 16:16:16, 48bpp, 16X, 16Y,
 355  *  16Z, with the 2-byte value for each X/Y/Z component stored as
 356  *  little-endian.
 357  */
 358 shared_ptr<dcp::OpenJPEGImage>
 359 dcp::xyz_to_xyz (uint8_t const * xyz_16, dcp::Size size, int stride)
 360 {
 361         shared_ptr<OpenJPEGImage> xyz_12 (new OpenJPEGImage (size));
 362
 363         int jn = 0;
 364         for (int y = 0; y < size.height; ++y) {
 365                 uint16_t const * p = reinterpret_cast<uint16_t const *> (xyz_16 + y * stride);
 366                 for (int x = 0; x < size.width; ++x) {
 367                         /* Truncate 16-bit to 12-bit */
 368                         xyz_12->data(0)[jn] = *p++ >> 4;
 369                         xyz_12->data(1)[jn] = *p++ >> 4;
 370                         xyz_12->data(2)[jn] = *p++ >> 4;
 371                         ++jn;
 372                 }
 373         }
 374
 375         return xyz_12;
 376 }