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 rgb RGB data; packed RGB 16:16:16, 48bpp, 16R, 16G, 16B,
 248  *  with the 2-byte value for each R/G/B component stored as
 249  *  little-endian; i.e. AV_PIX_FMT_RGB48LE.
 250  *  @param size of RGB image in pixels.
 251  *  @param stride of RGB data in pixels.
 252  */
 253 shared_ptr<dcp::OpenJPEGImage>
 254 dcp::rgb_to_xyz (
 255         uint8_t const * rgb,
 256         dcp::Size size,
 257         int stride,
 258         ColourConversion const & conversion,
 259         optional<NoteHandler> note
 260         )
 261 {
 262         shared_ptr<OpenJPEGImage> xyz (new OpenJPEGImage (size));
 263
 264         struct {
 265                 double r, g, b;
 266         } s;
 267
 268         struct {
 269                 double x, y, z;
 270         } d;
 271
 272         double const * lut_in = conversion.in()->lut (12, false);
 273         double const * lut_out = conversion.out()->lut (16, true);
 274         boost::numeric::ublas::matrix<double> const rgb_to_xyz = conversion.rgb_to_xyz ();
 275         boost::numeric::ublas::matrix<double> const bradford = conversion.bradford ();
 276
 277         /* This is is the product of the RGB to XYZ matrix, the Bradford transform and the DCI companding */
 278         double fast_matrix[9] = {
 279                 (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,
 280                 (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,
 281                 (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,
 282                 (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,
 283                 (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,
 284                 (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,
 285                 (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,
 286                 (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,
 287                 (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
 288                         };
 289
 290         int clamped = 0;
 291         int* xyz_x = xyz->data (0);
 292         int* xyz_y = xyz->data (1);
 293         int* xyz_z = xyz->data (2);
 294         for (int y = 0; y < size.height; ++y) {
 295                 uint16_t const * p = reinterpret_cast<uint16_t const *> (rgb + y * stride);
 296                 for (int x = 0; x < size.width; ++x) {
 297
 298                         /* In gamma LUT (converting 16-bit to 12-bit) */
 299                         s.r = lut_in[*p++ >> 4];
 300                         s.g = lut_in[*p++ >> 4];
 301                         s.b = lut_in[*p++ >> 4];
 302
 303                         /* RGB to XYZ, Bradford transform and DCI companding */
 304                         d.x = s.r * fast_matrix[0] + s.g * fast_matrix[1] + s.b * fast_matrix[2];
 305                         d.y = s.r * fast_matrix[3] + s.g * fast_matrix[4] + s.b * fast_matrix[5];
 306                         d.z = s.r * fast_matrix[6] + s.g * fast_matrix[7] + s.b * fast_matrix[8];
 307
 308                         /* Clamp */
 309
 310                         if (d.x < 0 || d.y < 0 || d.z < 0 || d.x > 65535 || d.y > 65535 || d.z > 65535) {
 311                                 ++clamped;
 312                         }
 313
 314                         d.x = max (0.0, d.x);
 315                         d.y = max (0.0, d.y);
 316                         d.z = max (0.0, d.z);
 317                         d.x = min (65535.0, d.x);
 318                         d.y = min (65535.0, d.y);
 319                         d.z = min (65535.0, d.z);
 320
 321                         /* Out gamma LUT */
 322                         *xyz_x++ = lrint (lut_out[lrint(d.x)] * 4095);
 323                         *xyz_y++ = lrint (lut_out[lrint(d.y)] * 4095);
 324                         *xyz_z++ = lrint (lut_out[lrint(d.z)] * 4095);
 325                 }
 326         }
 327
 328         if (clamped && note) {
 329                 note.get() (DCP_NOTE, String::compose ("%1 XYZ value(s) clamped", clamped));
 330         }
 331
 332         return xyz;
 333 }
 334
 335
 336 /** @param xyz_16 XYZ image data in packed 16:16:16, 48bpp, 16X, 16Y,
 337  *  16Z, with the 2-byte value for each X/Y/Z component stored as
 338  *  little-endian.
 339  */
 340 shared_ptr<dcp::OpenJPEGImage>
 341 dcp::xyz_to_xyz (uint8_t const * xyz_16, dcp::Size size, int stride)
 342 {
 343         shared_ptr<OpenJPEGImage> xyz_12 (new OpenJPEGImage (size));
 344
 345         int jn = 0;
 346         for (int y = 0; y < size.height; ++y) {
 347                 uint16_t const * p = reinterpret_cast<uint16_t const *> (xyz_16 + y * stride);
 348                 for (int x = 0; x < size.width; ++x) {
 349                         /* Truncate 16-bit to 12-bit */
 350                         xyz_12->data(0)[jn] = *p++ >> 4;
 351                         xyz_12->data(1)[jn] = *p++ >> 4;
 352                         xyz_12->data(2)[jn] = *p++ >> 4;
 353                         ++jn;
 354                 }
 355         }
 356
 357         return xyz_12;
 358 }