+ /* Otherwise one has an adjusted white and other hasn't, or they both have but different */
+ return false;
+}
+
+boost::numeric::ublas::matrix<double>
+ColourConversion::rgb_to_xyz () const
+{
+ /* See doc/design/colour.tex */
+
+ double const D = (_red.x - _white.x) * (_white.y - _blue.y) - (_white.x - _blue.x) * (_red.y - _white.y);
+ double const E = (_white.x - _green.x) * (_red.y - _white.y) - (_red.x - _white.x) * (_white.y - _green.y);
+ double const F = (_white.x - _green.x) * (_white.y - _blue.y) - (_white.x - _blue.x) * (_white.y - _green.y);
+ double const P = _red.y + _green.y * D / F + _blue.y * E / F;
+
+ boost::numeric::ublas::matrix<double> C (3, 3);
+ C(0, 0) = _red.x / P;
+ C(0, 1) = _green.x * D / (F * P);
+ C(0, 2) = _blue.x * E / (F * P);
+ C(1, 0) = _red.y / P;
+ C(1, 1) = _green.y * D / (F * P);
+ C(1, 2) = _blue.y * E / (F * P);
+ C(2, 0) = _red.z() / P;
+ C(2, 1) = _green.z() * D / (F * P);
+ C(2, 2) = _blue.z() * E / (F * P);
+ return C;
+}
+
+boost::numeric::ublas::matrix<double>
+ColourConversion::xyz_to_rgb () const
+{
+ boost::numeric::ublas::matrix<double> A (rgb_to_xyz ());
+
+ /* permutation matrix for the LU-factorization */
+ boost::numeric::ublas::permutation_matrix<std::size_t> pm (A.size1 ());
+
+ /* perform LU-factorization */
+ int const r = lu_factorize (A, pm);
+ DCP_ASSERT (r == 0);
+
+ /* create identity matrix of inverse */
+ boost::numeric::ublas::matrix<double> xyz_to_rgb (3, 3);
+ xyz_to_rgb.assign (boost::numeric::ublas::identity_matrix<double> (A.size1 ()));
+
+ /* backsubstitute to get the inverse */
+ lu_substitute (A, pm, xyz_to_rgb);
+
+ return xyz_to_rgb;
+}
+
+boost::numeric::ublas::matrix<double>
+ColourConversion::bradford () const
+{
+ if (!_adjusted_white || fabs (_adjusted_white.get().x) < 1e-6 || fabs (_adjusted_white.get().y) < 1e-6) {
+ boost::numeric::ublas::matrix<double> B = boost::numeric::ublas::zero_matrix<double> (3, 3);
+ B(0, 0) = 1;
+ B(1, 1) = 1;
+ B(2, 2) = 1;
+ return B;
+ }
+
+ /* See doc/design/colour.tex */
+
+ boost::numeric::ublas::matrix<double> M (3, 3);
+ M(0, 0) = 0.8951;
+ M(0, 1) = 0.2664;
+ M(0, 2) = -0.1614;
+ M(1, 0) = -0.7502;
+ M(1, 1) = 1.7135;
+ M(1, 2) = 0.0367;
+ M(2, 0) = 0.0389;
+ M(2, 1) = -0.0685;
+ M(2, 2) = 1.0296;
+
+ boost::numeric::ublas::matrix<double> Mi (3, 3);
+ Mi(0, 0) = 0.9869929055;
+ Mi(0, 1) = -0.1470542564;
+ Mi(0, 2) = 0.1599626517;
+ Mi(1, 0) = 0.4323052697;
+ Mi(1, 1) = 0.5183602715;
+ Mi(1, 2) = 0.0492912282;
+ Mi(2, 0) = -0.0085286646;
+ Mi(2, 1) = 0.0400428217;
+ Mi(2, 2) = 0.9684866958;
+
+ boost::numeric::ublas::matrix<double> Gp (3, 1);
+ Gp(0, 0) = _white.x / _white.y;
+ Gp(1, 0) = 1;
+ Gp(2, 0) = (1 - _white.x - _white.y) / _white.y;
+
+ boost::numeric::ublas::matrix<double> G = boost::numeric::ublas::prod (M, Gp);
+
+ boost::numeric::ublas::matrix<double> Hp (3, 1);
+ Hp(0, 0) = _adjusted_white.get().x / _adjusted_white.get().y;
+ Hp(1, 0) = 1;
+ Hp(2, 0) = (1 - _adjusted_white.get().x - _adjusted_white.get().y) / _adjusted_white.get().y;
+
+ boost::numeric::ublas::matrix<double> H = boost::numeric::ublas::prod (M, Hp);
+
+ boost::numeric::ublas::matrix<double> C = boost::numeric::ublas::zero_matrix<double> (3, 3);
+ C(0, 0) = H(0, 0) / G(0, 0);
+ C(1, 1) = H(1, 0) / G(1, 0);
+ C(2, 2) = H(2, 0) / G(2, 0);
+
+ boost::numeric::ublas::matrix<double> CM = boost::numeric::ublas::prod (C, M);
+ return boost::numeric::ublas::prod (Mi, CM);