src/colour_conversion.cc

   1 /*
   2     Copyright (C) 2014-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
  35 /** @file  src/colour_conversion.cc
  36  *  @brief ColourConversion class
  37  */
  38
  39
  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>
  49
  50
  51 using std::shared_ptr;
  52 using std::make_shared;
  53 using boost::optional;
  54 using namespace dcp;
  55
  56
  57 ColourConversion const &
  58 ColourConversion::srgb_to_xyz ()
  59 {
  60         static auto c = new ColourConversion (
  61                 make_shared<ModifiedGammaTransferFunction>(2.4, 0.04045, 0.055, 12.92),
  62                 YUVToRGB::REC601,
  63                 Chromaticity (0.64, 0.33),
  64                 Chromaticity (0.3, 0.6),
  65                 Chromaticity (0.15, 0.06),
  66                 Chromaticity::D65 (),
  67                 optional<Chromaticity> (),
  68                 make_shared<GammaTransferFunction>(2.6)
  69                 );
  70         return *c;
  71 }
  72
  73 ColourConversion const &
  74 ColourConversion::rec601_to_xyz ()
  75 {
  76         static auto c = new ColourConversion (
  77                 make_shared<GammaTransferFunction>(2.2),
  78                 YUVToRGB::REC601,
  79                 Chromaticity (0.64, 0.33),
  80                 Chromaticity (0.3, 0.6),
  81                 Chromaticity (0.15, 0.06),
  82                 Chromaticity::D65 (),
  83                 optional<Chromaticity> (),
  84                 make_shared<GammaTransferFunction>(2.6)
  85                 );
  86         return *c;
  87 }
  88
  89 ColourConversion const &
  90 ColourConversion::rec709_to_xyz ()
  91 {
  92         static auto c = new ColourConversion (
  93                 make_shared<GammaTransferFunction>(2.2),
  94                 YUVToRGB::REC709,
  95                 Chromaticity (0.64, 0.33),
  96                 Chromaticity (0.3, 0.6),
  97                 Chromaticity (0.15, 0.06),
  98                 Chromaticity::D65 (),
  99                 optional<Chromaticity> (),
 100                 make_shared<GammaTransferFunction>(2.6)
 101                 );
 102         return *c;
 103 }
 104
 105 ColourConversion const &
 106 ColourConversion::p3_to_xyz ()
 107 {
 108         static auto c = new ColourConversion (
 109                 make_shared<GammaTransferFunction>(2.6),
 110                 YUVToRGB::REC709,
 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)
 117                 );
 118         return *c;
 119 }
 120
 121 ColourConversion const &
 122 ColourConversion::rec1886_to_xyz ()
 123 {
 124         /* According to Olivier on DCP-o-matic bug #832, Rec. 1886 is Rec. 709 with
 125            2.4 gamma, so here goes ...
 126         */
 127         static auto c = new ColourConversion (
 128                 make_shared<GammaTransferFunction>(2.4),
 129                 YUVToRGB::REC709,
 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)
 136                 );
 137         return *c;
 138 }
 139
 140 ColourConversion const &
 141 ColourConversion::rec2020_to_xyz ()
 142 {
 143         static auto c = new ColourConversion (
 144                 make_shared<GammaTransferFunction>(2.4),
 145                 YUVToRGB::REC709,
 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)
 152                 );
 153         return *c;
 154 }
 155
 156 /** Sony S-Gamut3/S-Log3 */
 157 ColourConversion const &
 158 ColourConversion::s_gamut3_to_xyz ()
 159 {
 160         static auto c = new ColourConversion (
 161                 make_shared<SGamut3TransferFunction>(),
 162                 YUVToRGB::REC709,
 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>()
 169                 );
 170         return *c;
 171 }
 172
 173
 174
 175 ColourConversion::ColourConversion (
 176         shared_ptr<const TransferFunction> in,
 177         YUVToRGB yuv_to_rgb,
 178         Chromaticity red,
 179         Chromaticity green,
 180         Chromaticity blue,
 181         Chromaticity white,
 182         optional<Chromaticity> adjusted_white,
 183         shared_ptr<const TransferFunction> out
 184         )
 185         : _in (in)
 186         , _yuv_to_rgb (yuv_to_rgb)
 187         , _red (red)
 188         , _green (green)
 189         , _blue (blue)
 190         , _white (white)
 191         , _adjusted_white (adjusted_white)
 192         , _out (out)
 193 {
 194
 195 }
 196
 197
 198 bool
 199 ColourConversion::about_equal (ColourConversion const & other, float epsilon) const
 200 {
 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))) {
 210                 return false;
 211         }
 212
 213         if (!_adjusted_white && !other._adjusted_white) {
 214                 return true;
 215         }
 216
 217         if (
 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
 221                 ) {
 222                 return true;
 223         }
 224
 225         /* Otherwise one has an adjusted white and other hasn't, or they both have but different */
 226         return false;
 227 }
 228
 229
 230 boost::numeric::ublas::matrix<double>
 231 ColourConversion::rgb_to_xyz () const
 232 {
 233         /* See doc/design/colour.tex */
 234
 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;
 239
 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);
 250         return C;
 251 }
 252
 253
 254 boost::numeric::ublas::matrix<double>
 255 ColourConversion::xyz_to_rgb () const
 256 {
 257         boost::numeric::ublas::matrix<double> A (rgb_to_xyz());
 258
 259         /* permutation matrix for the LU-factorization */
 260         boost::numeric::ublas::permutation_matrix<std::size_t> pm (A.size1());
 261
 262         /* perform LU-factorization */
 263         int const r = lu_factorize (A, pm);
 264         DCP_ASSERT (r == 0);
 265
 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()));
 269
 270         /* backsubstitute to get the inverse */
 271         lu_substitute (A, pm, xyz_to_rgb);
 272
 273         return xyz_to_rgb;
 274 }
 275
 276
 277 boost::numeric::ublas::matrix<double>
 278 ColourConversion::bradford () const
 279 {
 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);
 282                 B(0, 0) = 1;
 283                 B(1, 1) = 1;
 284                 B(2, 2) = 1;
 285                 return B;
 286         }
 287
 288         /* See doc/design/colour.tex */
 289
 290         boost::numeric::ublas::matrix<double> M (3, 3);
 291         M(0, 0) = 0.8951;
 292         M(0, 1) = 0.2664;
 293         M(0, 2) = -0.1614;
 294         M(1, 0) = -0.7502;
 295         M(1, 1) = 1.7135;
 296         M(1, 2) = 0.0367;
 297         M(2, 0) = 0.0389;
 298         M(2, 1) = -0.0685;
 299         M(2, 2) = 1.0296;
 300
 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;
 311
 312         boost::numeric::ublas::matrix<double> Gp (3, 1);
 313         Gp(0, 0) = _white.x / _white.y;
 314         Gp(1, 0) = 1;
 315         Gp(2, 0) = (1 - _white.x - _white.y) / _white.y;
 316
 317         boost::numeric::ublas::matrix<double> G = boost::numeric::ublas::prod (M, Gp);
 318
 319         boost::numeric::ublas::matrix<double> Hp (3, 1);
 320         Hp(0, 0) = _adjusted_white.get().x / _adjusted_white.get().y;
 321         Hp(1, 0) = 1;
 322         Hp(2, 0) = (1 - _adjusted_white.get().x - _adjusted_white.get().y) / _adjusted_white.get().y;
 323
 324         boost::numeric::ublas::matrix<double> H = boost::numeric::ublas::prod (M, Hp);
 325
 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);
 330
 331         boost::numeric::ublas::matrix<double> CM = boost::numeric::ublas::prod (C, M);
 332         return boost::numeric::ublas::prod (Mi, CM);
 333 }