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