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 "colour_matrix.h"
  40 #include "dcp_assert.h"
  41 #include <boost/numeric/ublas/matrix.hpp>
  42 #include <boost/numeric/ublas/lu.hpp>
  43 #include <boost/numeric/ublas/io.hpp>
  44
  45 using boost::shared_ptr;
  46 using boost::optional;
  47 using namespace dcp;
  48
  49 ColourConversion const &
  50 ColourConversion::srgb_to_xyz ()
  51 {
  52         static ColourConversion* c = new ColourConversion (
  53                 shared_ptr<const TransferFunction> (new ModifiedGammaTransferFunction (2.4, 0.04045, 0.055, 12.92)),
  54                 YUV_TO_RGB_REC601,
  55                 Chromaticity (0.64, 0.33),
  56                 Chromaticity (0.3, 0.6),
  57                 Chromaticity (0.15, 0.06),
  58                 Chromaticity::D65 (),
  59                 optional<Chromaticity> (),
  60                 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6))
  61                 );
  62         return *c;
  63 }
  64
  65 ColourConversion const &
  66 ColourConversion::rec601_to_xyz ()
  67 {
  68         static ColourConversion* c = new ColourConversion (
  69                 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.2)),
  70                 YUV_TO_RGB_REC601,
  71                 Chromaticity (0.64, 0.33),
  72                 Chromaticity (0.3, 0.6),
  73                 Chromaticity (0.15, 0.06),
  74                 Chromaticity::D65 (),
  75                 optional<Chromaticity> (),
  76                 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6))
  77                 );
  78         return *c;
  79 }
  80
  81 ColourConversion const &
  82 ColourConversion::rec709_to_xyz ()
  83 {
  84         static ColourConversion* c = new ColourConversion (
  85                 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.2)),
  86                 YUV_TO_RGB_REC709,
  87                 Chromaticity (0.64, 0.33),
  88                 Chromaticity (0.3, 0.6),
  89                 Chromaticity (0.15, 0.06),
  90                 Chromaticity::D65 (),
  91                 optional<Chromaticity> (),
  92                 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6))
  93                 );
  94         return *c;
  95 }
  96
  97 ColourConversion const &
  98 ColourConversion::p3_to_xyz ()
  99 {
 100         static ColourConversion* c = new ColourConversion (
 101                 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6)),
 102                 YUV_TO_RGB_REC709,
 103                 Chromaticity (0.68, 0.32),
 104                 Chromaticity (0.265, 0.69),
 105                 Chromaticity (0.15, 0.06),
 106                 Chromaticity (0.314, 0.351),
 107                 optional<Chromaticity> (),
 108                 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6))
 109                 );
 110         return *c;
 111 }
 112
 113 ColourConversion const &
 114 ColourConversion::rec1886_to_xyz ()
 115 {
 116         /* According to Olivier on DCP-o-matic bug #832, Rec. 1886 is Rec. 709 with
 117            2.4 gamma, so here goes ...
 118         */
 119         static ColourConversion* c = new ColourConversion (
 120                 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.4)),
 121                 YUV_TO_RGB_REC709,
 122                 Chromaticity (0.64, 0.33),
 123                 Chromaticity (0.3, 0.6),
 124                 Chromaticity (0.15, 0.06),
 125                 Chromaticity::D65 (),
 126                 optional<Chromaticity> (),
 127                 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6))
 128                 );
 129         return *c;
 130 }
 131
 132 ColourConversion const &
 133 ColourConversion::rec2020_to_xyz ()
 134 {
 135         static ColourConversion* c = new ColourConversion (
 136                 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.4)),
 137                 YUV_TO_RGB_REC709,
 138                 Chromaticity (0.708, 0.292),
 139                 Chromaticity (0.170, 0.797),
 140                 Chromaticity (0.131, 0.046),
 141                 Chromaticity::D65 (),
 142                 optional<Chromaticity> (),
 143                 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6))
 144                 );
 145         return *c;
 146 }
 147
 148 /** Sony S-Gamut3/S-Log3 */
 149 ColourConversion const &
 150 ColourConversion::s_gamut3_to_xyz ()
 151 {
 152         static ColourConversion* c = new ColourConversion (
 153                 shared_ptr<const TransferFunction> (new SGamut3TransferFunction ()),
 154                 YUV_TO_RGB_REC709,
 155                 Chromaticity (0.73, 0.280),
 156                 Chromaticity (0.140, 0.855),
 157                 Chromaticity (0.100, -0.050),
 158                 Chromaticity::D65 (),
 159                 optional<Chromaticity> (),
 160                 shared_ptr<const TransferFunction> (new IdentityTransferFunction ())
 161                 );
 162         return *c;
 163 }
 164
 165
 166
 167 ColourConversion::ColourConversion (
 168         shared_ptr<const TransferFunction> in,
 169         YUVToRGB yuv_to_rgb,
 170         Chromaticity red,
 171         Chromaticity green,
 172         Chromaticity blue,
 173         Chromaticity white,
 174         optional<Chromaticity> adjusted_white,
 175         shared_ptr<const TransferFunction> out
 176         )
 177         : _in (in)
 178         , _yuv_to_rgb (yuv_to_rgb)
 179         , _red (red)
 180         , _green (green)
 181         , _blue (blue)
 182         , _white (white)
 183         , _adjusted_white (adjusted_white)
 184         , _out (out)
 185 {
 186
 187 }
 188
 189 bool
 190 ColourConversion::about_equal (ColourConversion const & other, float epsilon) const
 191 {
 192         if (!_in->about_equal (other._in, epsilon) ||
 193             _yuv_to_rgb != other._yuv_to_rgb ||
 194             !_red.about_equal (other._red, epsilon) ||
 195             !_green.about_equal (other._green, epsilon) ||
 196             !_blue.about_equal (other._blue, epsilon) ||
 197             !_white.about_equal (other._white, epsilon) ||
 198             (!_out && other._out) ||
 199             (_out && !other._out) ||
 200             (_out && !_out->about_equal (other._out, epsilon))) {
 201                 return false;
 202         }
 203
 204         if (!_adjusted_white && !other._adjusted_white) {
 205                 return true;
 206         }
 207
 208         if (
 209                 _adjusted_white && other._adjusted_white &&
 210                 fabs (_adjusted_white.get().x - other._adjusted_white.get().x) < epsilon &&
 211                 fabs (_adjusted_white.get().y - other._adjusted_white.get().y) < epsilon
 212                 ) {
 213                 return true;
 214         }
 215
 216         /* Otherwise one has an adjusted white and other hasn't, or they both have but different */
 217         return false;
 218 }
 219
 220 boost::numeric::ublas::matrix<double>
 221 ColourConversion::rgb_to_xyz () const
 222 {
 223         /* See doc/design/colour.tex */
 224
 225         double const D = (_red.x - _white.x) * (_white.y - _blue.y) - (_white.x - _blue.x) * (_red.y - _white.y);
 226         double const E = (_white.x - _green.x) * (_red.y - _white.y) - (_red.x - _white.x) * (_white.y - _green.y);
 227         double const F = (_white.x - _green.x) * (_white.y - _blue.y) - (_white.x - _blue.x) * (_white.y - _green.y);
 228         double const P = _red.y + _green.y * D / F + _blue.y * E / F;
 229
 230         boost::numeric::ublas::matrix<double> C (3, 3);
 231         C(0, 0) = _red.x / P;
 232         C(0, 1) = _green.x * D / (F * P);
 233         C(0, 2) = _blue.x * E / (F * P);
 234         C(1, 0) = _red.y / P;
 235         C(1, 1) = _green.y * D / (F * P);
 236         C(1, 2) = _blue.y * E / (F * P);
 237         C(2, 0) = _red.z() / P;
 238         C(2, 1) = _green.z() * D / (F * P);
 239         C(2, 2) = _blue.z() * E / (F * P);
 240         return C;
 241 }
 242
 243 boost::numeric::ublas::matrix<double>
 244 ColourConversion::xyz_to_rgb () const
 245 {
 246         boost::numeric::ublas::matrix<double> A (rgb_to_xyz ());
 247
 248         /* permutation matrix for the LU-factorization */
 249         boost::numeric::ublas::permutation_matrix<std::size_t> pm (A.size1 ());
 250
 251         /* perform LU-factorization */
 252         int const r = lu_factorize (A, pm);
 253         DCP_ASSERT (r == 0);
 254
 255         /* create identity matrix of inverse */
 256         boost::numeric::ublas::matrix<double> xyz_to_rgb (3, 3);
 257         xyz_to_rgb.assign (boost::numeric::ublas::identity_matrix<double> (A.size1 ()));
 258
 259         /* backsubstitute to get the inverse */
 260         lu_substitute (A, pm, xyz_to_rgb);
 261
 262         return xyz_to_rgb;
 263 }
 264
 265 boost::numeric::ublas::matrix<double>
 266 ColourConversion::bradford () const
 267 {
 268         if (!_adjusted_white || fabs (_adjusted_white.get().x) < 1e-6 || fabs (_adjusted_white.get().y) < 1e-6) {
 269                 boost::numeric::ublas::matrix<double> B = boost::numeric::ublas::zero_matrix<double> (3, 3);
 270                 B(0, 0) = 1;
 271                 B(1, 1) = 1;
 272                 B(2, 2) = 1;
 273                 return B;
 274         }
 275
 276         /* See doc/design/colour.tex */
 277
 278         boost::numeric::ublas::matrix<double> M (3, 3);
 279         M(0, 0) = 0.8951;
 280         M(0, 1) = 0.2664;
 281         M(0, 2) = -0.1614;
 282         M(1, 0) = -0.7502;
 283         M(1, 1) = 1.7135;
 284         M(1, 2) = 0.0367;
 285         M(2, 0) = 0.0389;
 286         M(2, 1) = -0.0685;
 287         M(2, 2) = 1.0296;
 288
 289         boost::numeric::ublas::matrix<double> Mi (3, 3);
 290         Mi(0, 0) = 0.9869929055;
 291         Mi(0, 1) = -0.1470542564;
 292         Mi(0, 2) = 0.1599626517;
 293         Mi(1, 0) = 0.4323052697;
 294         Mi(1, 1) = 0.5183602715;
 295         Mi(1, 2) = 0.0492912282;
 296         Mi(2, 0) = -0.0085286646;
 297         Mi(2, 1) = 0.0400428217;
 298         Mi(2, 2) = 0.9684866958;
 299
 300         boost::numeric::ublas::matrix<double> Gp (3, 1);
 301         Gp(0, 0) = _white.x / _white.y;
 302         Gp(1, 0) = 1;
 303         Gp(2, 0) = (1 - _white.x - _white.y) / _white.y;
 304
 305         boost::numeric::ublas::matrix<double> G = boost::numeric::ublas::prod (M, Gp);
 306
 307         boost::numeric::ublas::matrix<double> Hp (3, 1);
 308         Hp(0, 0) = _adjusted_white.get().x / _adjusted_white.get().y;
 309         Hp(1, 0) = 1;
 310         Hp(2, 0) = (1 - _adjusted_white.get().x - _adjusted_white.get().y) / _adjusted_white.get().y;
 311
 312         boost::numeric::ublas::matrix<double> H = boost::numeric::ublas::prod (M, Hp);
 313
 314         boost::numeric::ublas::matrix<double> C = boost::numeric::ublas::zero_matrix<double> (3, 3);
 315         C(0, 0) = H(0, 0) / G(0, 0);
 316         C(1, 1) = H(1, 0) / G(1, 0);
 317         C(2, 2) = H(2, 0) / G(2, 0);
 318
 319         boost::numeric::ublas::matrix<double> CM = boost::numeric::ublas::prod (C, M);
 320         return boost::numeric::ublas::prod (Mi, CM);
 321 }