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_conversion.h"
  37 #include "transfer_function.h"
  38 #include "dcp_assert.h"
  39 #include "compose.hpp"
  40 #include <cmath>
  41
  42 using std::min;
  43 using std::max;
  44 using std::cout;
  45 using std::shared_ptr;
  46 using boost::optional;
  47 using namespace dcp;
  48
  49 #define DCI_COEFFICIENT (48.0 / 52.37)
  50
  51 /** Convert an XYZ image to RGBA.
  52  *  @param xyz_image Image in XYZ.
  53  *  @param conversion Colour conversion to use.
  54  *  @param argb Buffer to fill with RGBA data.  The format of the data is:
  55  *
  56  *  <pre>
  57  *  Byte   /- 0 -------|- 1 --------|- 2 --------|- 3 --------|- 4 --------|- 5 --------| ...
  58  *         |(0, 0) Blue|(0, 0)Green |(0, 0) Red  |(0, 0) Alpha|(0, 1) Blue |(0, 1) Green| ...
  59  *  </pre>
  60  *
  61  *  So that the first byte is the blue component of the pixel at x=0, y=0, the second
  62  *  is the green component, and so on.
  63  *
  64  *  Lines are packed so that the second row directly follows the first.
  65  */
  66 void
  67 dcp::xyz_to_rgba (
  68         std::shared_ptr<const OpenJPEGImage> xyz_image,
  69         ColourConversion const & conversion,
  70         uint8_t* argb,
  71         int stride
  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                 argb += stride;
 139         }
 140 }
 141
 142 /** Convert an XYZ image to 48bpp RGB.
 143  *  @param xyz_image Frame in XYZ.
 144  *  @param conversion Colour conversion to use.
 145  *  @param rgb Buffer to fill with RGB data.  Format is packed RGB
 146  *  16:16:16, 48bpp, 16R, 16G, 16B, with the 2-byte value for each
 147  *  R/G/B component stored as little-endian; i.e. AV_PIX_FMT_RGB48LE.
 148  *  @param stride Stride for RGB data in bytes.
 149  *  @param note Optional handler for any notes that may be made during the conversion (e.g. when clamping occurs).
 150  */
 151 void
 152 dcp::xyz_to_rgb (
 153         shared_ptr<const OpenJPEGImage> xyz_image,
 154         ColourConversion const & conversion,
 155         uint8_t* rgb,
 156         int stride,
 157         optional<NoteHandler> note
 158         )
 159 {
 160         struct {
 161                 double x, y, z;
 162         } s;
 163
 164         struct {
 165                 double r, g, b;
 166         } d;
 167
 168         /* These should be 12-bit values from 0-4095 */
 169         int* xyz_x = xyz_image->data (0);
 170         int* xyz_y = xyz_image->data (1);
 171         int* xyz_z = xyz_image->data (2);
 172
 173         double const * lut_in = conversion.out()->lut (12, false);
 174         double const * lut_out = conversion.in()->lut (16, true);
 175         boost::numeric::ublas::matrix<double> const matrix = conversion.xyz_to_rgb ();
 176
 177         double fast_matrix[9] = {
 178                 matrix (0, 0), matrix (0, 1), matrix (0, 2),
 179                 matrix (1, 0), matrix (1, 1), matrix (1, 2),
 180                 matrix (2, 0), matrix (2, 1), matrix (2, 2)
 181         };
 182
 183         int const height = xyz_image->size().height;
 184         int const width = xyz_image->size().width;
 185
 186         for (int y = 0; y < height; ++y) {
 187                 uint16_t* rgb_line = reinterpret_cast<uint16_t*> (rgb + y * stride);
 188                 for (int x = 0; x < width; ++x) {
 189
 190                         int cx = *xyz_x++;
 191                         int cy = *xyz_y++;
 192                         int cz = *xyz_z++;
 193
 194                         if (cx < 0 || cx > 4095) {
 195                                 if (note) {
 196                                         note.get() (DCP_NOTE, String::compose ("XYZ value %1 out of range", cx));
 197                                 }
 198                                 cx = max (min (cx, 4095), 0);
 199                         }
 200
 201                         if (cy < 0 || cy > 4095) {
 202                                 if (note) {
 203                                         note.get() (DCP_NOTE, String::compose ("XYZ value %1 out of range", cy));
 204                                 }
 205                                 cy = max (min (cy, 4095), 0);
 206                         }
 207
 208                         if (cz < 0 || cz > 4095) {
 209                                 if (note) {
 210                                         note.get() (DCP_NOTE, String::compose ("XYZ value %1 out of range", cz));
 211                                 }
 212                                 cz = max (min (cz, 4095), 0);
 213                         }
 214
 215                         /* In gamma LUT */
 216                         s.x = lut_in[cx];
 217                         s.y = lut_in[cy];
 218                         s.z = lut_in[cz];
 219
 220                         /* DCI companding */
 221                         s.x /= DCI_COEFFICIENT;
 222                         s.y /= DCI_COEFFICIENT;
 223                         s.z /= DCI_COEFFICIENT;
 224
 225                         /* XYZ to RGB */
 226                         d.r = ((s.x * fast_matrix[0]) + (s.y * fast_matrix[1]) + (s.z * fast_matrix[2]));
 227                         d.g = ((s.x * fast_matrix[3]) + (s.y * fast_matrix[4]) + (s.z * fast_matrix[5]));
 228                         d.b = ((s.x * fast_matrix[6]) + (s.y * fast_matrix[7]) + (s.z * fast_matrix[8]));
 229
 230                         d.r = min (d.r, 1.0);
 231                         d.r = max (d.r, 0.0);
 232
 233                         d.g = min (d.g, 1.0);
 234                         d.g = max (d.g, 0.0);
 235
 236                         d.b = min (d.b, 1.0);
 237                         d.b = max (d.b, 0.0);
 238
 239                         *rgb_line++ = lrint(lut_out[lrint(d.r * 65535)] * 65535);
 240                         *rgb_line++ = lrint(lut_out[lrint(d.g * 65535)] * 65535);
 241                         *rgb_line++ = lrint(lut_out[lrint(d.b * 65535)] * 65535);
 242                 }
 243         }
 244 }
 245
 246 /** @param conversion Colour conversion.
 247  *  @param matrix Filled in with the product of the RGB to XYZ matrix, the Bradford transform and the DCI companding.
 248  */
 249 void
 250 dcp::combined_rgb_to_xyz (ColourConversion const & conversion, double* matrix)
 251 {
 252         boost::numeric::ublas::matrix<double> const rgb_to_xyz = conversion.rgb_to_xyz ();
 253         boost::numeric::ublas::matrix<double> const bradford = conversion.bradford ();
 254
 255         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))
 256                 * DCI_COEFFICIENT * 65535;
 257         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))
 258                 * DCI_COEFFICIENT * 65535;
 259         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))
 260                 * DCI_COEFFICIENT * 65535;
 261         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))
 262                 * DCI_COEFFICIENT * 65535;
 263         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))
 264                 * DCI_COEFFICIENT * 65535;
 265         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))
 266                 * DCI_COEFFICIENT * 65535;
 267         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))
 268                 * DCI_COEFFICIENT * 65535;
 269         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))
 270                 * DCI_COEFFICIENT * 65535;
 271         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))
 272                 * DCI_COEFFICIENT * 65535;
 273 }
 274
 275 /** @param rgb RGB data; packed RGB 16:16:16, 48bpp, 16R, 16G, 16B,
 276  *  with the 2-byte value for each R/G/B component stored as
 277  *  little-endian; i.e. AV_PIX_FMT_RGB48LE.
 278  *  @param size size of RGB image in pixels.
 279  *  @param size stride of RGB data in pixels.
 280  */
 281 shared_ptr<dcp::OpenJPEGImage>
 282 dcp::rgb_to_xyz (
 283         uint8_t const * rgb,
 284         dcp::Size size,
 285         int stride,
 286         ColourConversion const & conversion,
 287         optional<NoteHandler> note
 288         )
 289 {
 290         shared_ptr<OpenJPEGImage> xyz (new OpenJPEGImage (size));
 291
 292         struct {
 293                 double r, g, b;
 294         } s;
 295
 296         struct {
 297                 double x, y, z;
 298         } d;
 299
 300         double const * lut_in = conversion.in()->lut (12, false);
 301         double const * lut_out = conversion.out()->lut (16, true);
 302
 303         /* This is is the product of the RGB to XYZ matrix, the Bradford transform and the DCI companding */
 304         double fast_matrix[9];
 305         combined_rgb_to_xyz (conversion, fast_matrix);
 306
 307         int clamped = 0;
 308         int* xyz_x = xyz->data (0);
 309         int* xyz_y = xyz->data (1);
 310         int* xyz_z = xyz->data (2);
 311         for (int y = 0; y < size.height; ++y) {
 312                 uint16_t const * p = reinterpret_cast<uint16_t const *> (rgb + y * stride);
 313                 for (int x = 0; x < size.width; ++x) {
 314
 315                         /* In gamma LUT (converting 16-bit to 12-bit) */
 316                         s.r = lut_in[*p++ >> 4];
 317                         s.g = lut_in[*p++ >> 4];
 318                         s.b = lut_in[*p++ >> 4];
 319
 320                         /* RGB to XYZ, Bradford transform and DCI companding */
 321                         d.x = s.r * fast_matrix[0] + s.g * fast_matrix[1] + s.b * fast_matrix[2];
 322                         d.y = s.r * fast_matrix[3] + s.g * fast_matrix[4] + s.b * fast_matrix[5];
 323                         d.z = s.r * fast_matrix[6] + s.g * fast_matrix[7] + s.b * fast_matrix[8];
 324
 325                         /* Clamp */
 326
 327                         if (d.x < 0 || d.y < 0 || d.z < 0 || d.x > 65535 || d.y > 65535 || d.z > 65535) {
 328                                 ++clamped;
 329                         }
 330
 331                         d.x = max (0.0, d.x);
 332                         d.y = max (0.0, d.y);
 333                         d.z = max (0.0, d.z);
 334                         d.x = min (65535.0, d.x);
 335                         d.y = min (65535.0, d.y);
 336                         d.z = min (65535.0, d.z);
 337
 338                         /* Out gamma LUT */
 339                         *xyz_x++ = lrint (lut_out[lrint(d.x)] * 4095);
 340                         *xyz_y++ = lrint (lut_out[lrint(d.y)] * 4095);
 341                         *xyz_z++ = lrint (lut_out[lrint(d.z)] * 4095);
 342                 }
 343         }
 344
 345         if (clamped && note) {
 346                 note.get() (DCP_NOTE, String::compose ("%1 XYZ value(s) clamped", clamped));
 347         }
 348
 349         return xyz;
 350 }