src/colour_conversion.cc

   1 /*
   2     Copyright (C) 2014-2015 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 #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>
  43
  44 using std::shared_ptr;
  45 using boost::optional;
  46 using namespace dcp;
  47
  48 ColourConversion const &
  49 ColourConversion::srgb_to_xyz ()
  50 {
  51         static ColourConversion* c = new ColourConversion (
  52                 shared_ptr<const TransferFunction> (new ModifiedGammaTransferFunction (2.4, 0.04045, 0.055, 12.92)),
  53                 YUV_TO_RGB_REC601,
  54                 Chromaticity (0.64, 0.33),
  55                 Chromaticity (0.3, 0.6),
  56                 Chromaticity (0.15, 0.06),
  57                 Chromaticity::D65 (),
  58                 optional<Chromaticity> (),
  59                 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6))
  60                 );
  61         return *c;
  62 }
  63
  64 ColourConversion const &
  65 ColourConversion::rec601_to_xyz ()
  66 {
  67         static ColourConversion* c = new ColourConversion (
  68                 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.2)),
  69                 YUV_TO_RGB_REC601,
  70                 Chromaticity (0.64, 0.33),
  71                 Chromaticity (0.3, 0.6),
  72                 Chromaticity (0.15, 0.06),
  73                 Chromaticity::D65 (),
  74                 optional<Chromaticity> (),
  75                 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6))
  76                 );
  77         return *c;
  78 }
  79
  80 ColourConversion const &
  81 ColourConversion::rec709_to_xyz ()
  82 {
  83         static ColourConversion* c = new ColourConversion (
  84                 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.2)),
  85                 YUV_TO_RGB_REC709,
  86                 Chromaticity (0.64, 0.33),
  87                 Chromaticity (0.3, 0.6),
  88                 Chromaticity (0.15, 0.06),
  89                 Chromaticity::D65 (),
  90                 optional<Chromaticity> (),
  91                 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6))
  92                 );
  93         return *c;
  94 }
  95
  96 ColourConversion const &
  97 ColourConversion::p3_to_xyz ()
  98 {
  99         static ColourConversion* c = new ColourConversion (
 100                 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6)),
 101                 YUV_TO_RGB_REC709,
 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))
 108                 );
 109         return *c;
 110 }
 111
 112 ColourConversion const &
 113 ColourConversion::rec1886_to_xyz ()
 114 {
 115         /* According to Olivier on DCP-o-matic bug #832, Rec. 1886 is Rec. 709 with
 116            2.4 gamma, so here goes ...
 117         */
 118         static ColourConversion* c = new ColourConversion (
 119                 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.4)),
 120                 YUV_TO_RGB_REC709,
 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))
 127                 );
 128         return *c;
 129 }
 130
 131 ColourConversion const &
 132 ColourConversion::rec2020_to_xyz ()
 133 {
 134         static ColourConversion* c = new ColourConversion (
 135                 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.4)),
 136                 YUV_TO_RGB_REC709,
 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))
 143                 );
 144         return *c;
 145 }
 146
 147 /** Sony S-Gamut3/S-Log3 */
 148 ColourConversion const &
 149 ColourConversion::s_gamut3_to_xyz ()
 150 {
 151         static ColourConversion* c = new ColourConversion (
 152                 shared_ptr<const TransferFunction> (new SGamut3TransferFunction ()),
 153                 YUV_TO_RGB_REC709,
 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 ())
 160                 );
 161         return *c;
 162 }
 163
 164
 165
 166 ColourConversion::ColourConversion (
 167         shared_ptr<const TransferFunction> in,
 168         YUVToRGB yuv_to_rgb,
 169         Chromaticity red,
 170         Chromaticity green,
 171         Chromaticity blue,
 172         Chromaticity white,
 173         optional<Chromaticity> adjusted_white,
 174         shared_ptr<const TransferFunction> out
 175         )
 176         : _in (in)
 177         , _yuv_to_rgb (yuv_to_rgb)
 178         , _red (red)
 179         , _green (green)
 180         , _blue (blue)
 181         , _white (white)
 182         , _adjusted_white (adjusted_white)
 183         , _out (out)
 184 {
 185
 186 }
 187
 188 bool
 189 ColourConversion::about_equal (ColourConversion const & other, float epsilon) const
 190 {
 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))) {
 200                 return false;
 201         }
 202
 203         if (!_adjusted_white && !other._adjusted_white) {
 204                 return true;
 205         }
 206
 207         if (
 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
 211                 ) {
 212                 return true;
 213         }
 214
 215         /* Otherwise one has an adjusted white and other hasn't, or they both have but different */
 216         return false;
 217 }
 218
 219 boost::numeric::ublas::matrix<double>
 220 ColourConversion::rgb_to_xyz () const
 221 {
 222         /* See doc/design/colour.tex */
 223
 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;
 228
 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);
 239         return C;
 240 }
 241
 242 boost::numeric::ublas::matrix<double>
 243 ColourConversion::xyz_to_rgb () const
 244 {
 245         boost::numeric::ublas::matrix<double> A (rgb_to_xyz ());
 246
 247         /* permutation matrix for the LU-factorization */
 248         boost::numeric::ublas::permutation_matrix<std::size_t> pm (A.size1 ());
 249
 250         /* perform LU-factorization */
 251         int const r = lu_factorize (A, pm);
 252         DCP_ASSERT (r == 0);
 253
 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 ()));
 257
 258         /* backsubstitute to get the inverse */
 259         lu_substitute (A, pm, xyz_to_rgb);
 260
 261         return xyz_to_rgb;
 262 }
 263
 264 boost::numeric::ublas::matrix<double>
 265 ColourConversion::bradford () const
 266 {
 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);
 269                 B(0, 0) = 1;
 270                 B(1, 1) = 1;
 271                 B(2, 2) = 1;
 272                 return B;
 273         }
 274
 275         /* See doc/design/colour.tex */
 276
 277         boost::numeric::ublas::matrix<double> M (3, 3);
 278         M(0, 0) = 0.8951;
 279         M(0, 1) = 0.2664;
 280         M(0, 2) = -0.1614;
 281         M(1, 0) = -0.7502;
 282         M(1, 1) = 1.7135;
 283         M(1, 2) = 0.0367;
 284         M(2, 0) = 0.0389;
 285         M(2, 1) = -0.0685;
 286         M(2, 2) = 1.0296;
 287
 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;
 298
 299         boost::numeric::ublas::matrix<double> Gp (3, 1);
 300         Gp(0, 0) = _white.x / _white.y;
 301         Gp(1, 0) = 1;
 302         Gp(2, 0) = (1 - _white.x - _white.y) / _white.y;
 303
 304         boost::numeric::ublas::matrix<double> G = boost::numeric::ublas::prod (M, Gp);
 305
 306         boost::numeric::ublas::matrix<double> Hp (3, 1);
 307         Hp(0, 0) = _adjusted_white.get().x / _adjusted_white.get().y;
 308         Hp(1, 0) = 1;
 309         Hp(2, 0) = (1 - _adjusted_white.get().x - _adjusted_white.get().y) / _adjusted_white.get().y;
 310
 311         boost::numeric::ublas::matrix<double> H = boost::numeric::ublas::prod (M, Hp);
 312
 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);
 317
 318         boost::numeric::ublas::matrix<double> CM = boost::numeric::ublas::prod (C, M);
 319         return boost::numeric::ublas::prod (Mi, CM);
 320 }