+ // Recover z = 1 - x - y.
+ double z_R = 1.0 - x_R - y_R;
+ double z_G = 1.0 - x_G - y_G;
+ double z_B = 1.0 - x_B - y_B;
+
+ // We have, for each primary (example is with red):
+ //
+ // X_R / (X_R + Y_R + Z_R) = x_R
+ // Y_R / (X_R + Y_R + Z_R) = y_R
+ // Z_R / (X_R + Y_R + Z_R) = z_R
+ //
+ // Some algebraic fiddling yields (unsurprisingly):
+ //
+ // X_R = (x_R / y_R) Y_R
+ // Z_R = (z_R / y_R) Y_R
+ //
+ // We also know that since RGB=(1,1,1) should give us the
+ // D65 illuminant, we must have
+ //
+ // X_R + X_G + X_B = D65_X
+ // Y_R + Y_G + Y_B = D65_Y
+ // Z_R + Z_G + Z_B = D65_Z
+ //
+ // But since we already know how to express Y and Z by
+ // some constant multiple of X, this reduces to
+ //
+ // k1 Y_R + k2 Y_G + k3 Y_B = D65_X
+ // Y_R + Y_G + Y_B = D65_Y
+ // k4 Y_R + k5 Y_G + k6 Y_B = D65_Z
+ //
+ // Which we can solve for (Y_R, Y_G, Y_B) by inverting a 3x3 matrix.
+
+ Matrix3d temp;
+ temp(0,0) = x_R / y_R;
+ temp(0,1) = x_G / y_G;
+ temp(0,2) = x_B / y_B;
+
+ temp(1,0) = 1.0;
+ temp(1,1) = 1.0;
+ temp(1,2) = 1.0;
+
+ temp(2,0) = z_R / y_R;
+ temp(2,1) = z_G / y_G;
+ temp(2,2) = z_B / y_B;
+
+ Vector3d d65_XYZ(d65_X, d65_Y, d65_Z);
+ Vector3d Y_RGB = temp.inverse() * d65_XYZ;
+
+ // Now convert xyY -> XYZ.
+ double X_R = temp(0,0) * Y_RGB[0];
+ double Z_R = temp(2,0) * Y_RGB[0];
+
+ double X_G = temp(0,1) * Y_RGB[1];
+ double Z_G = temp(2,1) * Y_RGB[1];
+
+ double X_B = temp(0,2) * Y_RGB[2];
+ double Z_B = temp(2,2) * Y_RGB[2];
+
+ Matrix3d m;
+ m(0,0) = X_R; m(0,1) = X_G; m(0,2) = X_B;
+ m(1,0) = Y_RGB[0]; m(1,1) = Y_RGB[1]; m(1,2) = Y_RGB[2];
+ m(2,0) = Z_R; m(2,1) = Z_G; m(2,2) = Z_B;
+
+ return m;