libdcp
rgb_xyz.cc
Go to the documentation of this file.
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 
40 #include "colour_conversion.h"
41 #include "compose.hpp"
42 #include "dcp_assert.h"
43 #include "openjpeg_image.h"
44 #include "rgb_xyz.h"
45 #include "transfer_function.h"
46 #include <cmath>
47 
48 
49 using std::cout;
50 using std::make_shared;
51 using std::max;
52 using std::min;
53 using std::shared_ptr;
54 using boost::optional;
55 using namespace dcp;
56 
57 
58 static auto constexpr DCI_COEFFICIENT = 48.0 / 52.37;
59 
60 
61 void
63  std::shared_ptr<const OpenJPEGImage> xyz_image,
64  ColourConversion const & conversion,
65  uint8_t* argb,
66  int stride
67  )
68 {
69  int const max_colour = pow (2, 16) - 1;
70 
71  struct {
72  double x, y, z;
73  } s;
74 
75  struct {
76  double r, g, b;
77  } d;
78 
79  int* xyz_x = xyz_image->data (0);
80  int* xyz_y = xyz_image->data (1);
81  int* xyz_z = xyz_image->data (2);
82 
83  double const * lut_in = conversion.out()->lut (12, false);
84  double const * lut_out = conversion.in()->lut (16, true);
85  boost::numeric::ublas::matrix<double> const matrix = conversion.xyz_to_rgb ();
86 
87  double fast_matrix[9] = {
88  matrix (0, 0), matrix (0, 1), matrix (0, 2),
89  matrix (1, 0), matrix (1, 1), matrix (1, 2),
90  matrix (2, 0), matrix (2, 1), matrix (2, 2)
91  };
92 
93  int const height = xyz_image->size().height;
94  int const width = xyz_image->size().width;
95 
96  for (int y = 0; y < height; ++y) {
97  uint8_t* argb_line = argb;
98  for (int x = 0; x < width; ++x) {
99 
100  DCP_ASSERT (*xyz_x >= 0 && *xyz_y >= 0 && *xyz_z >= 0 && *xyz_x < 4096 && *xyz_y < 4096 && *xyz_z < 4096);
101 
102  /* In gamma LUT */
103  s.x = lut_in[*xyz_x++];
104  s.y = lut_in[*xyz_y++];
105  s.z = lut_in[*xyz_z++];
106 
107  /* DCI companding */
108  s.x /= DCI_COEFFICIENT;
109  s.y /= DCI_COEFFICIENT;
110  s.z /= DCI_COEFFICIENT;
111 
112  /* XYZ to RGB */
113  d.r = ((s.x * fast_matrix[0]) + (s.y * fast_matrix[1]) + (s.z * fast_matrix[2]));
114  d.g = ((s.x * fast_matrix[3]) + (s.y * fast_matrix[4]) + (s.z * fast_matrix[5]));
115  d.b = ((s.x * fast_matrix[6]) + (s.y * fast_matrix[7]) + (s.z * fast_matrix[8]));
116 
117  d.r = min (d.r, 1.0);
118  d.r = max (d.r, 0.0);
119 
120  d.g = min (d.g, 1.0);
121  d.g = max (d.g, 0.0);
122 
123  d.b = min (d.b, 1.0);
124  d.b = max (d.b, 0.0);
125 
126  /* Out gamma LUT */
127  *argb_line++ = lut_out[lrint(d.b * max_colour)] * 0xff;
128  *argb_line++ = lut_out[lrint(d.g * max_colour)] * 0xff;
129  *argb_line++ = lut_out[lrint(d.r * max_colour)] * 0xff;
130  *argb_line++ = 0xff;
131  }
132 
133  argb += stride;
134  }
135 }
136 
137 
138 void
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  auto 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  auto 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()(NoteType::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()(NoteType::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()(NoteType::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 void
234 dcp::combined_rgb_to_xyz (ColourConversion const & conversion, double* matrix)
235 {
236  auto const rgb_to_xyz = conversion.rgb_to_xyz ();
237  auto const bradford = conversion.bradford ();
238 
239  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))
240  * DCI_COEFFICIENT * 65535;
241  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))
242  * DCI_COEFFICIENT * 65535;
243  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))
244  * DCI_COEFFICIENT * 65535;
245  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))
246  * DCI_COEFFICIENT * 65535;
247  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))
248  * DCI_COEFFICIENT * 65535;
249  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))
250  * DCI_COEFFICIENT * 65535;
251  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))
252  * DCI_COEFFICIENT * 65535;
253  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))
254  * DCI_COEFFICIENT * 65535;
255  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))
256  * DCI_COEFFICIENT * 65535;
257 }
258 
259 
260 shared_ptr<dcp::OpenJPEGImage>
262  uint8_t const * rgb,
263  dcp::Size size,
264  int stride,
265  ColourConversion const & conversion,
266  optional<NoteHandler> note
267  )
268 {
269  auto xyz = make_shared<OpenJPEGImage>(size);
270 
271  struct {
272  double r, g, b;
273  } s;
274 
275  struct {
276  double x, y, z;
277  } d;
278 
279  auto const * lut_in = conversion.in()->lut (12, false);
280  auto const * lut_out = conversion.out()->lut (16, true);
281 
282  /* This is is the product of the RGB to XYZ matrix, the Bradford transform and the DCI companding */
283  double fast_matrix[9];
284  combined_rgb_to_xyz (conversion, fast_matrix);
285 
286  int clamped = 0;
287  int* xyz_x = xyz->data (0);
288  int* xyz_y = xyz->data (1);
289  int* xyz_z = xyz->data (2);
290  for (int y = 0; y < size.height; ++y) {
291  auto p = reinterpret_cast<uint16_t const *> (rgb + y * stride);
292  for (int x = 0; x < size.width; ++x) {
293 
294  /* In gamma LUT (converting 16-bit to 12-bit) */
295  s.r = lut_in[*p++ >> 4];
296  s.g = lut_in[*p++ >> 4];
297  s.b = lut_in[*p++ >> 4];
298 
299  /* RGB to XYZ, Bradford transform and DCI companding */
300  d.x = s.r * fast_matrix[0] + s.g * fast_matrix[1] + s.b * fast_matrix[2];
301  d.y = s.r * fast_matrix[3] + s.g * fast_matrix[4] + s.b * fast_matrix[5];
302  d.z = s.r * fast_matrix[6] + s.g * fast_matrix[7] + s.b * fast_matrix[8];
303 
304  /* Clamp */
305 
306  if (d.x < 0 || d.y < 0 || d.z < 0 || d.x > 65535 || d.y > 65535 || d.z > 65535) {
307  ++clamped;
308  }
309 
310  d.x = max (0.0, d.x);
311  d.y = max (0.0, d.y);
312  d.z = max (0.0, d.z);
313  d.x = min (65535.0, d.x);
314  d.y = min (65535.0, d.y);
315  d.z = min (65535.0, d.z);
316 
317  /* Out gamma LUT */
318  *xyz_x++ = lrint (lut_out[lrint(d.x)] * 4095);
319  *xyz_y++ = lrint (lut_out[lrint(d.y)] * 4095);
320  *xyz_z++ = lrint (lut_out[lrint(d.z)] * 4095);
321  }
322  }
323 
324  if (clamped && note) {
325  note.get()(NoteType::NOTE, String::compose("%1 XYZ value(s) clamped", clamped));
326  }
327 
328  return xyz;
329 }
A representation of all the parameters involved the colourspace conversion of a YUV image to XYZ (via...
ColourConversion class.
DCP_ASSERT macro.
Namespace for everything in libdcp.
Definition: array_data.h:50
void combined_rgb_to_xyz(ColourConversion const &conversion, double *matrix)
Definition: rgb_xyz.cc:234
void xyz_to_rgb(std::shared_ptr< const OpenJPEGImage >, ColourConversion const &conversion, uint8_t *rgb, int stride, boost::optional< NoteHandler > note=boost::optional< NoteHandler >())
std::shared_ptr< OpenJPEGImage > rgb_to_xyz(uint8_t const *rgb, dcp::Size size, int stride, ColourConversion const &conversion, boost::optional< NoteHandler > note=boost::optional< NoteHandler >())
void xyz_to_rgba(std::shared_ptr< const OpenJPEGImage >, ColourConversion const &conversion, uint8_t *rgba, int stride)
Definition: rgb_xyz.cc:62
OpenJPEGImage class.
Conversion between RGB and XYZ.
The integer, two-dimensional size of something.
Definition: types.h:71
TransferFunction class.