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         /* From Wikipedia */
 138         static ColourConversion* c = new ColourConversion (
 139                 shared_ptr<const TransferFunction> (new ModifiedGammaTransferFunction (1 / 0.45, 0.08145, 0.0993, 4.5)),
 140                 YUV_TO_RGB_REC709,
 141                 Chromaticity (0.708, 0.292),
 142                 Chromaticity (0.170, 0.797),
 143                 Chromaticity (0.131, 0.046),
 144                 /* D65 */
 145                 Chromaticity (0.3127, 0.329),
 146                 optional<Chromaticity> (),
 147                 shared_ptr<const TransferFunction> (new GammaTransferFunction (2.6))
 148                 );
 149         return *c;
 150 }
 151
 152 ColourConversion::ColourConversion (
 153         shared_ptr<const TransferFunction> in,
 154         YUVToRGB yuv_to_rgb,
 155         Chromaticity red,
 156         Chromaticity green,
 157         Chromaticity blue,
 158         Chromaticity white,
 159         optional<Chromaticity> adjusted_white,
 160         shared_ptr<const TransferFunction> out
 161         )
 162         : _in (in)
 163         , _yuv_to_rgb (yuv_to_rgb)
 164         , _red (red)
 165         , _green (green)
 166         , _blue (blue)
 167         , _white (white)
 168         , _adjusted_white (adjusted_white)
 169         , _out (out)
 170 {
 171
 172 }
 173
 174 bool
 175 ColourConversion::about_equal (ColourConversion const & other, float epsilon) const
 176 {
 177         if (!_in->about_equal (other._in, epsilon) ||
 178             _yuv_to_rgb != other._yuv_to_rgb ||
 179             !_red.about_equal (other._red, epsilon) ||
 180             !_green.about_equal (other._green, epsilon) ||
 181             !_blue.about_equal (other._blue, epsilon) ||
 182             !_white.about_equal (other._white, epsilon) ||
 183             !_out->about_equal (other._out, epsilon)) {
 184                 return false;
 185         }
 186
 187         if (!_adjusted_white && !other._adjusted_white) {
 188                 return true;
 189         }
 190
 191         if (
 192                 _adjusted_white && other._adjusted_white &&
 193                 fabs (_adjusted_white.get().x - other._adjusted_white.get().x) < epsilon &&
 194                 fabs (_adjusted_white.get().y - other._adjusted_white.get().y) < epsilon
 195                 ) {
 196                 return true;
 197         }
 198
 199         /* Otherwise one has an adjusted white and other hasn't, or they both have but different */
 200         return false;
 201 }
 202
 203 boost::numeric::ublas::matrix<double>
 204 ColourConversion::rgb_to_xyz () const
 205 {
 206         /* See doc/design/colour.tex */
 207
 208         double const D = (_red.x - _white.x) * (_white.y - _blue.y) - (_white.x - _blue.x) * (_red.y - _white.y);
 209         double const E = (_white.x - _green.x) * (_red.y - _white.y) - (_red.x - _white.x) * (_white.y - _green.y);
 210         double const F = (_white.x - _green.x) * (_white.y - _blue.y) - (_white.x - _blue.x) * (_white.y - _green.y);
 211         double const P = _red.y + _green.y * D / F + _blue.y * E / F;
 212
 213         boost::numeric::ublas::matrix<double> C (3, 3);
 214         C(0, 0) = _red.x / P;
 215         C(0, 1) = _green.x * D / (F * P);
 216         C(0, 2) = _blue.x * E / (F * P);
 217         C(1, 0) = _red.y / P;
 218         C(1, 1) = _green.y * D / (F * P);
 219         C(1, 2) = _blue.y * E / (F * P);
 220         C(2, 0) = _red.z() / P;
 221         C(2, 1) = _green.z() * D / (F * P);
 222         C(2, 2) = _blue.z() * E / (F * P);
 223         return C;
 224 }
 225
 226 boost::numeric::ublas::matrix<double>
 227 ColourConversion::xyz_to_rgb () const
 228 {
 229         boost::numeric::ublas::matrix<double> A (rgb_to_xyz ());
 230
 231         /* permutation matrix for the LU-factorization */
 232         boost::numeric::ublas::permutation_matrix<std::size_t> pm (A.size1 ());
 233
 234         /* perform LU-factorization */
 235         int const r = lu_factorize (A, pm);
 236         DCP_ASSERT (r == 0);
 237
 238         /* create identity matrix of inverse */
 239         boost::numeric::ublas::matrix<double> xyz_to_rgb (3, 3);
 240         xyz_to_rgb.assign (boost::numeric::ublas::identity_matrix<double> (A.size1 ()));
 241
 242         /* backsubstitute to get the inverse */
 243         lu_substitute (A, pm, xyz_to_rgb);
 244
 245         return xyz_to_rgb;
 246 }
 247
 248 boost::numeric::ublas::matrix<double>
 249 ColourConversion::bradford () const
 250 {
 251         if (!_adjusted_white || fabs (_adjusted_white.get().x) < 1e-6 || fabs (_adjusted_white.get().y) < 1e-6) {
 252                 boost::numeric::ublas::matrix<double> B = boost::numeric::ublas::zero_matrix<double> (3, 3);
 253                 B(0, 0) = 1;
 254                 B(1, 1) = 1;
 255                 B(2, 2) = 1;
 256                 return B;
 257         }
 258
 259         /* See doc/design/colour.tex */
 260
 261         boost::numeric::ublas::matrix<double> M (3, 3);
 262         M(0, 0) = 0.8951;
 263         M(0, 1) = 0.2664;
 264         M(0, 2) = -0.1614;
 265         M(1, 0) = -0.7502;
 266         M(1, 1) = 1.7135;
 267         M(1, 2) = 0.0367;
 268         M(2, 0) = 0.0389;
 269         M(2, 1) = -0.0685;
 270         M(2, 2) = 1.0296;
 271
 272         boost::numeric::ublas::matrix<double> Mi (3, 3);
 273         Mi(0, 0) = 0.9869929055;
 274         Mi(0, 1) = -0.1470542564;
 275         Mi(0, 2) = 0.1599626517;
 276         Mi(1, 0) = 0.4323052697;
 277         Mi(1, 1) = 0.5183602715;
 278         Mi(1, 2) = 0.0492912282;
 279         Mi(2, 0) = -0.0085286646;
 280         Mi(2, 1) = 0.0400428217;
 281         Mi(2, 2) = 0.9684866958;
 282
 283         boost::numeric::ublas::matrix<double> Gp (3, 1);
 284         Gp(0, 0) = _white.x / _white.y;
 285         Gp(1, 0) = 1;
 286         Gp(2, 0) = (1 - _white.x - _white.y) / _white.y;
 287
 288         boost::numeric::ublas::matrix<double> G = boost::numeric::ublas::prod (M, Gp);
 289
 290         boost::numeric::ublas::matrix<double> Hp (3, 1);
 291         Hp(0, 0) = _adjusted_white.get().x / _adjusted_white.get().y;
 292         Hp(1, 0) = 1;
 293         Hp(2, 0) = (1 - _adjusted_white.get().x - _adjusted_white.get().y) / _adjusted_white.get().y;
 294
 295         boost::numeric::ublas::matrix<double> H = boost::numeric::ublas::prod (M, Hp);
 296
 297         boost::numeric::ublas::matrix<double> C = boost::numeric::ublas::zero_matrix<double> (3, 3);
 298         C(0, 0) = H(0, 0) / G(0, 0);
 299         C(1, 1) = H(1, 0) / G(1, 0);
 300         C(2, 2) = H(2, 0) / G(2, 0);
 301
 302         boost::numeric::ublas::matrix<double> CM = boost::numeric::ublas::prod (C, M);
 303         return boost::numeric::ublas::prod (Mi, CM);
 304 }