2 Copyright (C) 2014-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.
35 /** @file src/colour_conversion.cc
36 * @brief ColourConversion class
40 #include "colour_conversion.h"
41 #include "gamma_transfer_function.h"
42 #include "modified_gamma_transfer_function.h"
43 #include "s_gamut3_transfer_function.h"
44 #include "identity_transfer_function.h"
45 #include "dcp_assert.h"
46 #include <boost/numeric/ublas/matrix.hpp>
47 #include <boost/numeric/ublas/lu.hpp>
48 #include <boost/numeric/ublas/io.hpp>
51 using std::shared_ptr;
52 using std::make_shared;
53 using boost::optional;
57 ColourConversion const &
58 ColourConversion::srgb_to_xyz ()
60 static auto c = new ColourConversion (
61 make_shared<ModifiedGammaTransferFunction>(2.4, 0.04045, 0.055, 12.92),
63 Chromaticity (0.64, 0.33),
64 Chromaticity (0.3, 0.6),
65 Chromaticity (0.15, 0.06),
67 optional<Chromaticity> (),
68 make_shared<GammaTransferFunction>(2.6)
73 ColourConversion const &
74 ColourConversion::rec601_to_xyz ()
76 static auto c = new ColourConversion (
77 make_shared<GammaTransferFunction>(2.2),
79 Chromaticity (0.64, 0.33),
80 Chromaticity (0.3, 0.6),
81 Chromaticity (0.15, 0.06),
83 optional<Chromaticity> (),
84 make_shared<GammaTransferFunction>(2.6)
89 ColourConversion const &
90 ColourConversion::rec709_to_xyz ()
92 static auto c = new ColourConversion (
93 make_shared<GammaTransferFunction>(2.2),
95 Chromaticity (0.64, 0.33),
96 Chromaticity (0.3, 0.6),
97 Chromaticity (0.15, 0.06),
99 optional<Chromaticity> (),
100 make_shared<GammaTransferFunction>(2.6)
105 ColourConversion const &
106 ColourConversion::p3_to_xyz ()
108 static auto c = new ColourConversion (
109 make_shared<GammaTransferFunction>(2.6),
111 Chromaticity (0.68, 0.32),
112 Chromaticity (0.265, 0.69),
113 Chromaticity (0.15, 0.06),
114 Chromaticity (0.314, 0.351),
115 optional<Chromaticity> (),
116 make_shared<GammaTransferFunction>(2.6)
121 ColourConversion const &
122 ColourConversion::rec1886_to_xyz ()
124 /* According to Olivier on DCP-o-matic bug #832, Rec. 1886 is Rec. 709 with
125 2.4 gamma, so here goes ...
127 static auto c = new ColourConversion (
128 make_shared<GammaTransferFunction>(2.4),
130 Chromaticity (0.64, 0.33),
131 Chromaticity (0.3, 0.6),
132 Chromaticity (0.15, 0.06),
133 Chromaticity::D65 (),
134 optional<Chromaticity> (),
135 make_shared<GammaTransferFunction>(2.6)
140 ColourConversion const &
141 ColourConversion::rec2020_to_xyz ()
143 static auto c = new ColourConversion (
144 make_shared<GammaTransferFunction>(2.4),
146 Chromaticity (0.708, 0.292),
147 Chromaticity (0.170, 0.797),
148 Chromaticity (0.131, 0.046),
149 Chromaticity::D65 (),
150 optional<Chromaticity> (),
151 make_shared<GammaTransferFunction>(2.6)
156 /** Sony S-Gamut3/S-Log3 */
157 ColourConversion const &
158 ColourConversion::s_gamut3_to_xyz ()
160 static auto c = new ColourConversion (
161 make_shared<SGamut3TransferFunction>(),
163 Chromaticity (0.73, 0.280),
164 Chromaticity (0.140, 0.855),
165 Chromaticity (0.100, -0.050),
166 Chromaticity::D65 (),
167 optional<Chromaticity> (),
168 make_shared<IdentityTransferFunction>()
175 ColourConversion::ColourConversion (
176 shared_ptr<const TransferFunction> in,
182 optional<Chromaticity> adjusted_white,
183 shared_ptr<const TransferFunction> out
186 , _yuv_to_rgb (yuv_to_rgb)
191 , _adjusted_white (adjusted_white)
199 ColourConversion::about_equal (ColourConversion const & other, float epsilon) const
201 if (!_in->about_equal (other._in, epsilon) ||
202 _yuv_to_rgb != other._yuv_to_rgb ||
203 !_red.about_equal (other._red, epsilon) ||
204 !_green.about_equal (other._green, epsilon) ||
205 !_blue.about_equal (other._blue, epsilon) ||
206 !_white.about_equal (other._white, epsilon) ||
207 (!_out && other._out) ||
208 (_out && !other._out) ||
209 (_out && !_out->about_equal (other._out, epsilon))) {
213 if (!_adjusted_white && !other._adjusted_white) {
218 _adjusted_white && other._adjusted_white &&
219 fabs (_adjusted_white.get().x - other._adjusted_white.get().x) < epsilon &&
220 fabs (_adjusted_white.get().y - other._adjusted_white.get().y) < epsilon
225 /* Otherwise one has an adjusted white and other hasn't, or they both have but different */
230 boost::numeric::ublas::matrix<double>
231 ColourConversion::rgb_to_xyz () const
233 /* See doc/design/colour.tex */
235 double const D = (_red.x - _white.x) * (_white.y - _blue.y) - (_white.x - _blue.x) * (_red.y - _white.y);
236 double const E = (_white.x - _green.x) * (_red.y - _white.y) - (_red.x - _white.x) * (_white.y - _green.y);
237 double const F = (_white.x - _green.x) * (_white.y - _blue.y) - (_white.x - _blue.x) * (_white.y - _green.y);
238 double const P = _red.y + _green.y * D / F + _blue.y * E / F;
240 boost::numeric::ublas::matrix<double> C (3, 3);
241 C(0, 0) = _red.x / P;
242 C(0, 1) = _green.x * D / (F * P);
243 C(0, 2) = _blue.x * E / (F * P);
244 C(1, 0) = _red.y / P;
245 C(1, 1) = _green.y * D / (F * P);
246 C(1, 2) = _blue.y * E / (F * P);
247 C(2, 0) = _red.z() / P;
248 C(2, 1) = _green.z() * D / (F * P);
249 C(2, 2) = _blue.z() * E / (F * P);
254 boost::numeric::ublas::matrix<double>
255 ColourConversion::xyz_to_rgb () const
257 boost::numeric::ublas::matrix<double> A (rgb_to_xyz());
259 /* permutation matrix for the LU-factorization */
260 boost::numeric::ublas::permutation_matrix<std::size_t> pm (A.size1());
262 /* perform LU-factorization */
263 int const r = lu_factorize (A, pm);
266 /* create identity matrix of inverse */
267 boost::numeric::ublas::matrix<double> xyz_to_rgb (3, 3);
268 xyz_to_rgb.assign (boost::numeric::ublas::identity_matrix<double> (A.size1()));
270 /* backsubstitute to get the inverse */
271 lu_substitute (A, pm, xyz_to_rgb);
277 boost::numeric::ublas::matrix<double>
278 ColourConversion::bradford () const
280 if (!_adjusted_white || fabs (_adjusted_white.get().x) < 1e-6 || fabs (_adjusted_white.get().y) < 1e-6) {
281 boost::numeric::ublas::matrix<double> B = boost::numeric::ublas::zero_matrix<double> (3, 3);
288 /* See doc/design/colour.tex */
290 boost::numeric::ublas::matrix<double> M (3, 3);
301 boost::numeric::ublas::matrix<double> Mi (3, 3);
302 Mi(0, 0) = 0.9869929055;
303 Mi(0, 1) = -0.1470542564;
304 Mi(0, 2) = 0.1599626517;
305 Mi(1, 0) = 0.4323052697;
306 Mi(1, 1) = 0.5183602715;
307 Mi(1, 2) = 0.0492912282;
308 Mi(2, 0) = -0.0085286646;
309 Mi(2, 1) = 0.0400428217;
310 Mi(2, 2) = 0.9684866958;
312 boost::numeric::ublas::matrix<double> Gp (3, 1);
313 Gp(0, 0) = _white.x / _white.y;
315 Gp(2, 0) = (1 - _white.x - _white.y) / _white.y;
317 boost::numeric::ublas::matrix<double> G = boost::numeric::ublas::prod (M, Gp);
319 boost::numeric::ublas::matrix<double> Hp (3, 1);
320 Hp(0, 0) = _adjusted_white.get().x / _adjusted_white.get().y;
322 Hp(2, 0) = (1 - _adjusted_white.get().x - _adjusted_white.get().y) / _adjusted_white.get().y;
324 boost::numeric::ublas::matrix<double> H = boost::numeric::ublas::prod (M, Hp);
326 boost::numeric::ublas::matrix<double> C = boost::numeric::ublas::zero_matrix<double> (3, 3);
327 C(0, 0) = H(0, 0) / G(0, 0);
328 C(1, 1) = H(1, 0) / G(1, 0);
329 C(2, 2) = H(2, 0) / G(2, 0);
331 boost::numeric::ublas::matrix<double> CM = boost::numeric::ublas::prod (C, M);
332 return boost::numeric::ublas::prod (Mi, CM);