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"
44 #include "piecewise_lut.h"
46 #include "transfer_function.h"
51 using std::make_shared;
54 using std::shared_ptr;
55 using boost::optional;
59 static auto constexpr DCI_COEFFICIENT = 48.0 / 52.37;
64 std::shared_ptr<const OpenJPEGImage> xyz_image,
65 ColourConversion const & conversion,
70 int const max_colour = pow (2, 16) - 1;
80 int* xyz_x = xyz_image->data (0);
81 int* xyz_y = xyz_image->data (1);
82 int* xyz_z = xyz_image->data (2);
84 auto lut_in = conversion.out()->double_lut(0, 1, 12, false);
85 auto lut_out = conversion.in()->double_lut(0, 1, 16, true);
86 boost::numeric::ublas::matrix<double> const matrix = conversion.xyz_to_rgb ();
88 double fast_matrix[9] = {
89 matrix (0, 0), matrix (0, 1), matrix (0, 2),
90 matrix (1, 0), matrix (1, 1), matrix (1, 2),
91 matrix (2, 0), matrix (2, 1), matrix (2, 2)
94 int const height = xyz_image->size().height;
95 int const width = xyz_image->size().width;
97 for (int y = 0; y < height; ++y) {
98 uint8_t* argb_line = argb;
99 for (int x = 0; x < width; ++x) {
101 DCP_ASSERT (*xyz_x >= 0 && *xyz_y >= 0 && *xyz_z >= 0 && *xyz_x < 4096 && *xyz_y < 4096 && *xyz_z < 4096);
104 s.x = lut_in[*xyz_x++];
105 s.y = lut_in[*xyz_y++];
106 s.z = lut_in[*xyz_z++];
109 s.x /= DCI_COEFFICIENT;
110 s.y /= DCI_COEFFICIENT;
111 s.z /= DCI_COEFFICIENT;
114 d.r = ((s.x * fast_matrix[0]) + (s.y * fast_matrix[1]) + (s.z * fast_matrix[2]));
115 d.g = ((s.x * fast_matrix[3]) + (s.y * fast_matrix[4]) + (s.z * fast_matrix[5]));
116 d.b = ((s.x * fast_matrix[6]) + (s.y * fast_matrix[7]) + (s.z * fast_matrix[8]));
118 d.r = min (d.r, 1.0);
119 d.r = max (d.r, 0.0);
121 d.g = min (d.g, 1.0);
122 d.g = max (d.g, 0.0);
124 d.b = min (d.b, 1.0);
125 d.b = max (d.b, 0.0);
128 *argb_line++ = lut_out[lrint(d.b * max_colour)] * 0xff;
129 *argb_line++ = lut_out[lrint(d.g * max_colour)] * 0xff;
130 *argb_line++ = lut_out[lrint(d.r * max_colour)] * 0xff;
141 shared_ptr<const OpenJPEGImage> xyz_image,
142 ColourConversion const & conversion,
145 optional<NoteHandler> note
156 /* These should be 12-bit values from 0-4095 */
157 int* xyz_x = xyz_image->data (0);
158 int* xyz_y = xyz_image->data (1);
159 int* xyz_z = xyz_image->data (2);
161 auto lut_in = conversion.out()->double_lut(0, 1, 12, false);
162 auto lut_out = conversion.in()->double_lut(0, 1, 16, true);
163 auto const matrix = conversion.xyz_to_rgb ();
165 double fast_matrix[9] = {
166 matrix (0, 0), matrix (0, 1), matrix (0, 2),
167 matrix (1, 0), matrix (1, 1), matrix (1, 2),
168 matrix (2, 0), matrix (2, 1), matrix (2, 2)
171 int const height = xyz_image->size().height;
172 int const width = xyz_image->size().width;
174 for (int y = 0; y < height; ++y) {
175 auto rgb_line = reinterpret_cast<uint16_t*> (rgb + y * stride);
176 for (int x = 0; x < width; ++x) {
182 if (cx < 0 || cx > 4095) {
184 note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cx));
186 cx = max (min (cx, 4095), 0);
189 if (cy < 0 || cy > 4095) {
191 note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cy));
193 cy = max (min (cy, 4095), 0);
196 if (cz < 0 || cz > 4095) {
198 note.get()(NoteType::NOTE, String::compose("XYZ value %1 out of range", cz));
200 cz = max (min (cz, 4095), 0);
209 s.x /= DCI_COEFFICIENT;
210 s.y /= DCI_COEFFICIENT;
211 s.z /= DCI_COEFFICIENT;
214 d.r = ((s.x * fast_matrix[0]) + (s.y * fast_matrix[1]) + (s.z * fast_matrix[2]));
215 d.g = ((s.x * fast_matrix[3]) + (s.y * fast_matrix[4]) + (s.z * fast_matrix[5]));
216 d.b = ((s.x * fast_matrix[6]) + (s.y * fast_matrix[7]) + (s.z * fast_matrix[8]));
218 d.r = min (d.r, 1.0);
219 d.r = max (d.r, 0.0);
221 d.g = min (d.g, 1.0);
222 d.g = max (d.g, 0.0);
224 d.b = min (d.b, 1.0);
225 d.b = max (d.b, 0.0);
227 *rgb_line++ = lrint(lut_out[lrint(d.r * 65535)] * 65535);
228 *rgb_line++ = lrint(lut_out[lrint(d.g * 65535)] * 65535);
229 *rgb_line++ = lrint(lut_out[lrint(d.b * 65535)] * 65535);
235 dcp::combined_rgb_to_xyz (ColourConversion const & conversion, double* matrix)
237 auto const rgb_to_xyz = conversion.rgb_to_xyz ();
238 auto const bradford = conversion.bradford ();
240 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))
242 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))
244 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))
246 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))
248 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))
250 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))
252 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))
254 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))
256 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))
262 dcp::make_inverse_gamma_lut(shared_ptr<const TransferFunction> fn)
264 /* The parameters here were chosen by trial and error to reduce errors when running rgb_xyz_lut_test */
265 return PiecewiseLUT2(fn, 0.062, 16, 12, true, 4095);
278 ColourConversion const& conversion
289 auto lut_in = conversion.in()->double_lut(0, 1, 12, false);
290 auto lut_out = make_inverse_gamma_lut(conversion.out());
292 /* This is is the product of the RGB to XYZ matrix, the Bradford transform and the DCI companding */
293 double fast_matrix[9];
294 combined_rgb_to_xyz (conversion, fast_matrix);
296 for (int y = 0; y < size.height; ++y) {
297 auto p = reinterpret_cast<uint16_t const *> (rgb + y * stride);
298 for (int x = 0; x < size.width; ++x) {
300 /* In gamma LUT (converting 16-bit to 12-bit) */
301 s.r = lut_in[*p++ >> 4];
302 s.g = lut_in[*p++ >> 4];
303 s.b = lut_in[*p++ >> 4];
305 /* RGB to XYZ, Bradford transform and DCI companding */
306 d.x = s.r * fast_matrix[0] + s.g * fast_matrix[1] + s.b * fast_matrix[2];
307 d.y = s.r * fast_matrix[3] + s.g * fast_matrix[4] + s.b * fast_matrix[5];
308 d.z = s.r * fast_matrix[6] + s.g * fast_matrix[7] + s.b * fast_matrix[8];
311 d.x = max (0.0, d.x);
312 d.y = max (0.0, d.y);
313 d.z = max (0.0, d.z);
314 d.x = min (1.0, d.x);
315 d.y = min (1.0, d.y);
316 d.z = min (1.0, d.z);
319 *xyz_x++ = lut_out.lookup(d.x);
320 *xyz_y++ = lut_out.lookup(d.y);
321 *xyz_z++ = lut_out.lookup(d.z);
327 shared_ptr<dcp::OpenJPEGImage>
332 ColourConversion const & conversion
335 auto xyz = make_shared<OpenJPEGImage>(size);
337 int* xyz_x = xyz->data (0);
338 int* xyz_y = xyz->data (1);
339 int* xyz_z = xyz->data (2);
341 rgb_to_xyz_internal(rgb, xyz_x, xyz_y, xyz_z, size, stride, conversion);
353 ColourConversion const & conversion
356 rgb_to_xyz_internal(rgb, dst, dst, dst, size, stride, conversion);