2 Copyright (C) 2014-2015 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.
34 #include "colour_conversion.h"
35 #include "gamma_transfer_function.h"
36 #include "modified_gamma_transfer_function.h"
37 #include "s_gamut3_transfer_function.h"
38 #include "identity_transfer_function.h"
39 #include "dcp_assert.h"
40 #include <boost/numeric/ublas/matrix.hpp>
41 #include <boost/numeric/ublas/lu.hpp>
42 #include <boost/numeric/ublas/io.hpp>
44 using boost::shared_ptr;
45 using boost::optional;
48 ColourConversion const &
49 ColourConversion::srgb_to_xyz ()
51 static ColourConversion* c = new ColourConversion (
52 shared_ptr<const TransferFunction> (new ModifiedGammaTransferFunction (2.4, 0.04045, 0.055, 12.92)),
54 Chromaticity (0.64, 0.33),
55 Chromaticity (0.3, 0.6),
56 Chromaticity (0.15, 0.06),
58 optional<Chromaticity> (),
59 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6))
64 ColourConversion const &
65 ColourConversion::rec601_to_xyz ()
67 static ColourConversion* c = new ColourConversion (
68 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.2)),
70 Chromaticity (0.64, 0.33),
71 Chromaticity (0.3, 0.6),
72 Chromaticity (0.15, 0.06),
74 optional<Chromaticity> (),
75 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6))
80 ColourConversion const &
81 ColourConversion::rec709_to_xyz ()
83 static ColourConversion* c = new ColourConversion (
84 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.2)),
86 Chromaticity (0.64, 0.33),
87 Chromaticity (0.3, 0.6),
88 Chromaticity (0.15, 0.06),
90 optional<Chromaticity> (),
91 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6))
96 ColourConversion const &
97 ColourConversion::p3_to_xyz ()
99 static ColourConversion* c = new ColourConversion (
100 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6)),
102 Chromaticity (0.68, 0.32),
103 Chromaticity (0.265, 0.69),
104 Chromaticity (0.15, 0.06),
105 Chromaticity (0.314, 0.351),
106 optional<Chromaticity> (),
107 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6))
112 ColourConversion const &
113 ColourConversion::rec1886_to_xyz ()
115 /* According to Olivier on DCP-o-matic bug #832, Rec. 1886 is Rec. 709 with
116 2.4 gamma, so here goes ...
118 static ColourConversion* c = new ColourConversion (
119 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.4)),
121 Chromaticity (0.64, 0.33),
122 Chromaticity (0.3, 0.6),
123 Chromaticity (0.15, 0.06),
124 Chromaticity::D65 (),
125 optional<Chromaticity> (),
126 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6))
131 ColourConversion const &
132 ColourConversion::rec2020_to_xyz ()
134 static ColourConversion* c = new ColourConversion (
135 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.4)),
137 Chromaticity (0.708, 0.292),
138 Chromaticity (0.170, 0.797),
139 Chromaticity (0.131, 0.046),
140 Chromaticity::D65 (),
141 optional<Chromaticity> (),
142 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6))
147 /** Sony S-Gamut3/S-Log3 */
148 ColourConversion const &
149 ColourConversion::s_gamut3_to_xyz ()
151 static ColourConversion* c = new ColourConversion (
152 shared_ptr<const TransferFunction> (new SGamut3TransferFunction ()),
154 Chromaticity (0.73, 0.280),
155 Chromaticity (0.140, 0.855),
156 Chromaticity (0.100, -0.050),
157 Chromaticity::D65 (),
158 optional<Chromaticity> (),
159 shared_ptr<const TransferFunction> (new IdentityTransferFunction ())
166 ColourConversion::ColourConversion (
167 shared_ptr<const TransferFunction> in,
173 optional<Chromaticity> adjusted_white,
174 shared_ptr<const TransferFunction> out
177 , _yuv_to_rgb (yuv_to_rgb)
182 , _adjusted_white (adjusted_white)
189 ColourConversion::about_equal (ColourConversion const & other, float epsilon) const
191 if (!_in->about_equal (other._in, epsilon) ||
192 _yuv_to_rgb != other._yuv_to_rgb ||
193 !_red.about_equal (other._red, epsilon) ||
194 !_green.about_equal (other._green, epsilon) ||
195 !_blue.about_equal (other._blue, epsilon) ||
196 !_white.about_equal (other._white, epsilon) ||
197 (!_out && other._out) ||
198 (_out && !other._out) ||
199 (_out && !_out->about_equal (other._out, epsilon))) {
203 if (!_adjusted_white && !other._adjusted_white) {
208 _adjusted_white && other._adjusted_white &&
209 fabs (_adjusted_white.get().x - other._adjusted_white.get().x) < epsilon &&
210 fabs (_adjusted_white.get().y - other._adjusted_white.get().y) < epsilon
215 /* Otherwise one has an adjusted white and other hasn't, or they both have but different */
219 boost::numeric::ublas::matrix<double>
220 ColourConversion::rgb_to_xyz () const
222 /* See doc/design/colour.tex */
224 double const D = (_red.x - _white.x) * (_white.y - _blue.y) - (_white.x - _blue.x) * (_red.y - _white.y);
225 double const E = (_white.x - _green.x) * (_red.y - _white.y) - (_red.x - _white.x) * (_white.y - _green.y);
226 double const F = (_white.x - _green.x) * (_white.y - _blue.y) - (_white.x - _blue.x) * (_white.y - _green.y);
227 double const P = _red.y + _green.y * D / F + _blue.y * E / F;
229 boost::numeric::ublas::matrix<double> C (3, 3);
230 C(0, 0) = _red.x / P;
231 C(0, 1) = _green.x * D / (F * P);
232 C(0, 2) = _blue.x * E / (F * P);
233 C(1, 0) = _red.y / P;
234 C(1, 1) = _green.y * D / (F * P);
235 C(1, 2) = _blue.y * E / (F * P);
236 C(2, 0) = _red.z() / P;
237 C(2, 1) = _green.z() * D / (F * P);
238 C(2, 2) = _blue.z() * E / (F * P);
242 boost::numeric::ublas::matrix<double>
243 ColourConversion::xyz_to_rgb () const
245 boost::numeric::ublas::matrix<double> A (rgb_to_xyz ());
247 /* permutation matrix for the LU-factorization */
248 boost::numeric::ublas::permutation_matrix<std::size_t> pm (A.size1 ());
250 /* perform LU-factorization */
251 int const r = lu_factorize (A, pm);
254 /* create identity matrix of inverse */
255 boost::numeric::ublas::matrix<double> xyz_to_rgb (3, 3);
256 xyz_to_rgb.assign (boost::numeric::ublas::identity_matrix<double> (A.size1 ()));
258 /* backsubstitute to get the inverse */
259 lu_substitute (A, pm, xyz_to_rgb);
264 boost::numeric::ublas::matrix<double>
265 ColourConversion::bradford () const
267 if (!_adjusted_white || fabs (_adjusted_white.get().x) < 1e-6 || fabs (_adjusted_white.get().y) < 1e-6) {
268 boost::numeric::ublas::matrix<double> B = boost::numeric::ublas::zero_matrix<double> (3, 3);
275 /* See doc/design/colour.tex */
277 boost::numeric::ublas::matrix<double> M (3, 3);
288 boost::numeric::ublas::matrix<double> Mi (3, 3);
289 Mi(0, 0) = 0.9869929055;
290 Mi(0, 1) = -0.1470542564;
291 Mi(0, 2) = 0.1599626517;
292 Mi(1, 0) = 0.4323052697;
293 Mi(1, 1) = 0.5183602715;
294 Mi(1, 2) = 0.0492912282;
295 Mi(2, 0) = -0.0085286646;
296 Mi(2, 1) = 0.0400428217;
297 Mi(2, 2) = 0.9684866958;
299 boost::numeric::ublas::matrix<double> Gp (3, 1);
300 Gp(0, 0) = _white.x / _white.y;
302 Gp(2, 0) = (1 - _white.x - _white.y) / _white.y;
304 boost::numeric::ublas::matrix<double> G = boost::numeric::ublas::prod (M, Gp);
306 boost::numeric::ublas::matrix<double> Hp (3, 1);
307 Hp(0, 0) = _adjusted_white.get().x / _adjusted_white.get().y;
309 Hp(2, 0) = (1 - _adjusted_white.get().x - _adjusted_white.get().y) / _adjusted_white.get().y;
311 boost::numeric::ublas::matrix<double> H = boost::numeric::ublas::prod (M, Hp);
313 boost::numeric::ublas::matrix<double> C = boost::numeric::ublas::zero_matrix<double> (3, 3);
314 C(0, 0) = H(0, 0) / G(0, 0);
315 C(1, 1) = H(1, 0) / G(1, 0);
316 C(2, 2) = H(2, 0) / G(2, 0);
318 boost::numeric::ublas::matrix<double> CM = boost::numeric::ublas::prod (C, M);
319 return boost::numeric::ublas::prod (Mi, CM);