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