src/rgb_xyz.cc

   1 /*
   2     Copyright (C) 2013-2015 Carl Hetherington <cth@carlh.net>
   3
   4     This program is free software; you can redistribute it and/or modify
   5     it under the terms of the GNU General Public License as published by
   6     the Free Software Foundation; either version 2 of the License, or
   7     (at your option) any later version.
   8
   9     This program is distributed in the hope that it will be useful,
  10     but WITHOUT ANY WARRANTY; without even the implied warranty of
  11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12     GNU General Public License for more details.
  13
  14     You should have received a copy of the GNU General Public License
  15     along with this program; if not, write to the Free Software
  16     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  17
  18 */
  19
  20 #include "rgb_xyz.h"
  21 #include "openjpeg_image.h"
  22 #include "colour_matrix.h"
  23 #include "colour_conversion.h"
  24 #include "transfer_function.h"
  25 #include "dcp_assert.h"
  26 #include "compose.hpp"
  27 #include <cmath>
  28
  29 using std::min;
  30 using std::max;
  31 using std::cout;
  32 using boost::shared_ptr;
  33 using boost::optional;
  34 using namespace dcp;
  35
  36 #define DCI_COEFFICIENT (48.0 / 52.37)
  37
  38 /** Convert an XYZ image to RGBA.
  39  *  @param xyz_image Image in XYZ.
  40  *  @param conversion Colour conversion to use.
  41  *  @param argb Buffer to fill with RGBA data.  The format of the data is:
  42  *
  43  *  <pre>
  44  *  Byte   /- 0 -------|- 1 --------|- 2 --------|- 3 --------|- 4 --------|- 5 --------| ...
  45  *         |(0, 0) Blue|(0, 0)Green |(0, 0) Red  |(0, 0) Alpha|(0, 1) Blue |(0, 1) Green| ...
  46  *  </pre>
  47  *
  48  *  So that the first byte is the blue component of the pixel at x=0, y=0, the second
  49  *  is the green component, and so on.
  50  *
  51  *  Lines are packed so that the second row directly follows the first.
  52  */
  53 void
  54 dcp::xyz_to_rgba (
  55         boost::shared_ptr<const OpenJPEGImage> xyz_image,
  56         ColourConversion const & conversion,
  57         uint8_t* argb
  58         )
  59 {
  60         int const max_colour = pow (2, 12) - 1;
  61
  62         struct {
  63                 double x, y, z;
  64         } s;
  65
  66         struct {
  67                 double r, g, b;
  68         } d;
  69
  70         int* xyz_x = xyz_image->data (0);
  71         int* xyz_y = xyz_image->data (1);
  72         int* xyz_z = xyz_image->data (2);
  73
  74         double const * lut_in = conversion.out()->lut (12, false);
  75         double const * lut_out = conversion.in()->lut (16, true);
  76         boost::numeric::ublas::matrix<double> const matrix = conversion.xyz_to_rgb ();
  77
  78         double fast_matrix[9] = {
  79                 matrix (0, 0), matrix (0, 1), matrix (0, 2),
  80                 matrix (1, 0), matrix (1, 1), matrix (1, 2),
  81                 matrix (2, 0), matrix (2, 1), matrix (2, 2)
  82         };
  83
  84         int const height = xyz_image->size().height;
  85         int const width = xyz_image->size().width;
  86
  87         for (int y = 0; y < height; ++y) {
  88                 uint8_t* argb_line = argb;
  89                 for (int x = 0; x < width; ++x) {
  90
  91                         DCP_ASSERT (*xyz_x >= 0 && *xyz_y >= 0 && *xyz_z >= 0 && *xyz_x < 4096 && *xyz_y < 4096 && *xyz_z < 4096);
  92
  93                         /* In gamma LUT */
  94                         s.x = lut_in[*xyz_x++];
  95                         s.y = lut_in[*xyz_y++];
  96                         s.z = lut_in[*xyz_z++];
  97
  98                         /* DCI companding */
  99                         s.x /= DCI_COEFFICIENT;
 100                         s.y /= DCI_COEFFICIENT;
 101                         s.z /= DCI_COEFFICIENT;
 102
 103                         /* XYZ to RGB */
 104                         d.r = ((s.x * fast_matrix[0]) + (s.y * fast_matrix[1]) + (s.z * fast_matrix[2]));
 105                         d.g = ((s.x * fast_matrix[3]) + (s.y * fast_matrix[4]) + (s.z * fast_matrix[5]));
 106                         d.b = ((s.x * fast_matrix[6]) + (s.y * fast_matrix[7]) + (s.z * fast_matrix[8]));
 107
 108                         d.r = min (d.r, 1.0);
 109                         d.r = max (d.r, 0.0);
 110
 111                         d.g = min (d.g, 1.0);
 112                         d.g = max (d.g, 0.0);
 113
 114                         d.b = min (d.b, 1.0);
 115                         d.b = max (d.b, 0.0);
 116
 117                         /* Out gamma LUT */
 118                         *argb_line++ = lut_out[lrint(d.b * max_colour)] * 0xff;
 119                         *argb_line++ = lut_out[lrint(d.g * max_colour)] * 0xff;
 120                         *argb_line++ = lut_out[lrint(d.r * max_colour)] * 0xff;
 121                         *argb_line++ = 0xff;
 122                 }
 123
 124                 /* 4 bytes per pixel */
 125                 argb += width * 4;
 126         }
 127 }
 128
 129 /** Convert an XYZ image to 48bpp RGB.
 130  *  @param xyz_image Frame in XYZ.
 131  *  @param conversion Colour conversion to use.
 132  *  @param rgb Buffer to fill with RGB data.  Format is packed RGB
 133  *  16:16:16, 48bpp, 16R, 16G, 16B, with the 2-byte value for each
 134  *  R/G/B component stored as little-endian; i.e. AV_PIX_FMT_RGB48LE.
 135  *  @param stride Stride for RGB data in bytes.
 136  *  @param note Optional handler for any notes that may be made during the conversion (e.g. when clamping occurs).
 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         boost::numeric::ublas::matrix<double> 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                 uint16_t* 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() (DCP_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() (DCP_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() (DCP_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 /** @param rgb RGB data; packed RGB 16:16:16, 48bpp, 16R, 16G, 16B,
 234  *  with the 2-byte value for each R/G/B component stored as
 235  *  little-endian; i.e. AV_PIX_FMT_RGB48LE.
 236  *  @param size of RGB image in pixels.
 237  *  @param stride of RGB data in pixels.
 238  */
 239 shared_ptr<dcp::OpenJPEGImage>
 240 dcp::rgb_to_xyz (
 241         uint8_t const * rgb,
 242         dcp::Size size,
 243         int stride,
 244         ColourConversion const & conversion,
 245         optional<NoteHandler> note
 246         )
 247 {
 248         shared_ptr<OpenJPEGImage> xyz (new OpenJPEGImage (size));
 249
 250         struct {
 251                 double r, g, b;
 252         } s;
 253
 254         struct {
 255                 double x, y, z;
 256         } d;
 257
 258         double const * lut_in = conversion.in()->lut (12, false);
 259         double const * lut_out = conversion.out()->lut (16, true);
 260         boost::numeric::ublas::matrix<double> const rgb_to_xyz = conversion.rgb_to_xyz ();
 261         boost::numeric::ublas::matrix<double> const bradford = conversion.bradford ();
 262
 263         /* This is is the product of the RGB to XYZ matrix, the Bradford transform and the DCI companding */
 264         double fast_matrix[9] = {
 265                 (bradford (0, 0) * rgb_to_xyz (0, 0) + bradford (0, 1) * rgb_to_xyz (1, 0) + bradford (0, 2) * rgb_to_xyz (2, 0)) * DCI_COEFFICIENT * 65535,
 266                 (bradford (0, 0) * rgb_to_xyz (0, 1) + bradford (0, 1) * rgb_to_xyz (1, 1) + bradford (0, 2) * rgb_to_xyz (2, 1)) * DCI_COEFFICIENT * 65535,
 267                 (bradford (0, 0) * rgb_to_xyz (0, 2) + bradford (0, 1) * rgb_to_xyz (1, 2) + bradford (0, 2) * rgb_to_xyz (2, 2)) * DCI_COEFFICIENT * 65535,
 268                 (bradford (1, 0) * rgb_to_xyz (0, 0) + bradford (1, 1) * rgb_to_xyz (1, 0) + bradford (1, 2) * rgb_to_xyz (2, 0)) * DCI_COEFFICIENT * 65535,
 269                 (bradford (1, 0) * rgb_to_xyz (0, 1) + bradford (1, 1) * rgb_to_xyz (1, 1) + bradford (1, 2) * rgb_to_xyz (2, 1)) * DCI_COEFFICIENT * 65535,
 270                 (bradford (1, 0) * rgb_to_xyz (0, 2) + bradford (1, 1) * rgb_to_xyz (1, 2) + bradford (1, 2) * rgb_to_xyz (2, 2)) * DCI_COEFFICIENT * 65535,
 271                 (bradford (2, 0) * rgb_to_xyz (0, 0) + bradford (2, 1) * rgb_to_xyz (1, 0) + bradford (2, 2) * rgb_to_xyz (2, 0)) * DCI_COEFFICIENT * 65535,
 272                 (bradford (2, 0) * rgb_to_xyz (0, 1) + bradford (2, 1) * rgb_to_xyz (1, 1) + bradford (2, 2) * rgb_to_xyz (2, 1)) * DCI_COEFFICIENT * 65535,
 273                 (bradford (2, 0) * rgb_to_xyz (0, 2) + bradford (2, 1) * rgb_to_xyz (1, 2) + bradford (2, 2) * rgb_to_xyz (2, 2)) * DCI_COEFFICIENT * 65535
 274                         };
 275
 276         int clamped = 0;
 277         int* xyz_x = xyz->data (0);
 278         int* xyz_y = xyz->data (1);
 279         int* xyz_z = xyz->data (2);
 280         for (int y = 0; y < size.height; ++y) {
 281                 uint16_t const * p = reinterpret_cast<uint16_t const *> (rgb + y * stride);
 282                 for (int x = 0; x < size.width; ++x) {
 283
 284                         /* In gamma LUT (converting 16-bit to 12-bit) */
 285                         s.r = lut_in[*p++ >> 4];
 286                         s.g = lut_in[*p++ >> 4];
 287                         s.b = lut_in[*p++ >> 4];
 288
 289                         /* RGB to XYZ, Bradford transform and DCI companding */
 290                         d.x = s.r * fast_matrix[0] + s.g * fast_matrix[1] + s.b * fast_matrix[2];
 291                         d.y = s.r * fast_matrix[3] + s.g * fast_matrix[4] + s.b * fast_matrix[5];
 292                         d.z = s.r * fast_matrix[6] + s.g * fast_matrix[7] + s.b * fast_matrix[8];
 293
 294                         /* Clamp */
 295
 296                         if (d.x < 0 || d.y < 0 || d.z < 0 || d.x > 65535 || d.y > 65535 || d.z > 65535) {
 297                                 ++clamped;
 298                         }
 299
 300                         d.x = max (0.0, d.x);
 301                         d.y = max (0.0, d.y);
 302                         d.z = max (0.0, d.z);
 303                         d.x = min (65535.0, d.x);
 304                         d.y = min (65535.0, d.y);
 305                         d.z = min (65535.0, d.z);
 306
 307                         /* Out gamma LUT */
 308                         *xyz_x++ = lut_out[lrint(d.x)] * 4095;
 309                         *xyz_y++ = lut_out[lrint(d.y)] * 4095;
 310                         *xyz_z++ = lut_out[lrint(d.z)] * 4095;
 311                 }
 312         }
 313
 314         if (clamped && note) {
 315                 note.get() (DCP_NOTE, String::compose ("%1 XYZ value(s) clamped", clamped));
 316         }
 317
 318         return xyz;
 319 }
 320
 321
 322 /** @param xyz_16 XYZ image data in packed 16:16:16, 48bpp, 16X, 16Y,
 323  *  16Z, with the 2-byte value for each X/Y/Z component stored as
 324  *  little-endian.
 325  */
 326 shared_ptr<dcp::OpenJPEGImage>
 327 dcp::xyz_to_xyz (uint8_t const * xyz_16, dcp::Size size, int stride)
 328 {
 329         shared_ptr<OpenJPEGImage> xyz_12 (new OpenJPEGImage (size));
 330
 331         int jn = 0;
 332         for (int y = 0; y < size.height; ++y) {
 333                 uint16_t const * p = reinterpret_cast<uint16_t const *> (xyz_16 + y * stride);
 334                 for (int x = 0; x < size.width; ++x) {
 335                         /* Truncate 16-bit to 12-bit */
 336                         xyz_12->data(0)[jn] = *p++ >> 4;
 337                         xyz_12->data(1)[jn] = *p++ >> 4;
 338                         xyz_12->data(2)[jn] = *p++ >> 4;
 339                         ++jn;
 340                 }
 341         }
 342
 343         return xyz_12;
 344 }