2 Copyright (C) 2013-2021 Carl Hetherington <cth@carlh.net>
4 This file is part of libdcp.
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.
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.
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/>.
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
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.
36 * @brief Conversion between RGB and XYZ
40 #include "colour_conversion.h"
41 #include "compose.hpp"
42 #include "dcp_assert.h"
43 #include "openjpeg_image.h"
45 #include "transfer_function.h"
50 using std::make_shared;
53 using std::shared_ptr;
54 using boost::optional;
58 static auto constexpr DCI_COEFFICIENT = 48.0 / 52.37;
63 std::shared_ptr<const OpenJPEGImage> xyz_image,
64 ColourConversion const & conversion,
69 int const max_colour = pow (2, 16) - 1;
79 int* xyz_x = xyz_image->data (0);
80 int* xyz_y = xyz_image->data (1);
81 int* xyz_z = xyz_image->data (2);
83 auto lut_in = conversion.out()->lut (12, false);
84 auto lut_out = conversion.in()->lut (16, true);
85 boost::numeric::ublas::matrix<double> const matrix = conversion.xyz_to_rgb ();
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)
93 int const height = xyz_image->size().height;
94 int const width = xyz_image->size().width;
96 for (int y = 0; y < height; ++y) {
97 uint8_t* argb_line = argb;
98 for (int x = 0; x < width; ++x) {
100 DCP_ASSERT (*xyz_x >= 0 && *xyz_y >= 0 && *xyz_z >= 0 && *xyz_x < 4096 && *xyz_y < 4096 && *xyz_z < 4096);
103 s.x = lut_in[*xyz_x++];
104 s.y = lut_in[*xyz_y++];
105 s.z = lut_in[*xyz_z++];
108 s.x /= DCI_COEFFICIENT;
109 s.y /= DCI_COEFFICIENT;
110 s.z /= DCI_COEFFICIENT;
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]));
117 d.r = min (d.r, 1.0);
118 d.r = max (d.r, 0.0);
120 d.g = min (d.g, 1.0);
121 d.g = max (d.g, 0.0);
123 d.b = min (d.b, 1.0);
124 d.b = max (d.b, 0.0);
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;
140 shared_ptr<const OpenJPEGImage> xyz_image,
141 ColourConversion const & conversion,
144 optional<NoteHandler> note
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);
160 auto lut_in = conversion.out()->lut (12, false);
161 auto lut_out = conversion.in()->lut (16, true);
162 auto const matrix = conversion.xyz_to_rgb ();
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)
170 int const height = xyz_image->size().height;
171 int const width = xyz_image->size().width;
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) {
181 if (cx < 0 || cx > 4095) {
183 note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cx));
185 cx = max (min (cx, 4095), 0);
188 if (cy < 0 || cy > 4095) {
190 note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cy));
192 cy = max (min (cy, 4095), 0);
195 if (cz < 0 || cz > 4095) {
197 note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cz));
199 cz = max (min (cz, 4095), 0);
208 s.x /= DCI_COEFFICIENT;
209 s.y /= DCI_COEFFICIENT;
210 s.z /= DCI_COEFFICIENT;
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]));
217 d.r = min (d.r, 1.0);
218 d.r = max (d.r, 0.0);
220 d.g = min (d.g, 1.0);
221 d.g = max (d.g, 0.0);
223 d.b = min (d.b, 1.0);
224 d.b = max (d.b, 0.0);
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);
234 dcp::combined_rgb_to_xyz (ColourConversion const & conversion, double* matrix)
236 auto const rgb_to_xyz = conversion.rgb_to_xyz ();
237 auto const bradford = conversion.bradford ();
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;
260 shared_ptr<dcp::OpenJPEGImage>
265 ColourConversion const & conversion,
266 optional<NoteHandler> note
269 auto xyz = make_shared<OpenJPEGImage>(size);
279 auto lut_in = conversion.in()->lut (12, false);
280 auto lut_out = conversion.out()->lut (16, true);
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);
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) {
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];
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];
306 if (d.x < 0 || d.y < 0 || d.z < 0 || d.x > 65535 || d.y > 65535 || d.z > 65535) {
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);
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);
324 if (clamped && note) {
325 note.get()(NoteType::NOTE, String::compose("%1 XYZ value(s) clamped", clamped));