+ dest[i] = src[i] / info->slope;
+}
+
+static inline void
+apply_matrix(double dest[3], double src[3], double matrix[9])
+{
+ dest[0] = src[0] * matrix[0] + src[1] * matrix[1] + src[2] * matrix[2];
+ dest[1] = src[0] * matrix[3] + src[1] * matrix[4] + src[2] * matrix[5];
+ dest[2] = src[0] * matrix[6] + src[1] * matrix[7] + src[2] * matrix[8];
+}
+
+void
+color_invert_matrix(double dest[9], double matrix[9])
+{
+ double *i = dest, *m = matrix;
+ double a0 = m[4] * m[8] - m[5] * m[7];
+ double a1 = m[3] * m[8] - m[5] * m[6];
+ double a2 = m[3] * m[7] - m[4] * m[6];
+ double d = 1 / (m[0] * a0 - m[1] * a1 + m[2] * a2);
+ i[0] = d * a0;
+ i[3] = -d * a1;
+ i[6] = d * a2;
+ i[1] = -d * (m[1] * m[8] - m[2] * m[7]);
+ i[4] = d * (m[0] * m[8] - m[2] * m[6]);
+ i[7] = -d * (m[0] * m[7] - m[1] * m[6]);
+ i[2] = d * (m[1] * m[5] - m[2] * m[4]);
+ i[5] = -d * (m[0] * m[5] - m[2] * m[3]);
+ i[8] = d * (m[0] * m[4] - m[1] * m[3]);
+}
+
+static void
+mul_matrices(double r[9], double a[9], double b[9])
+{
+ r[0] = a[0] * b[0] + a[1] * b[3] + a[2] * b[6];
+ r[1] = a[0] * b[1] + a[1] * b[4] + a[2] * b[7];
+ r[2] = a[0] * b[2] + a[1] * b[5] + a[2] * b[8];
+ r[3] = a[3] * b[0] + a[4] * b[3] + a[5] * b[6];
+ r[4] = a[3] * b[1] + a[4] * b[4] + a[5] * b[7];
+ r[5] = a[3] * b[2] + a[4] * b[5] + a[5] * b[8];
+ r[6] = a[6] * b[0] + a[7] * b[3] + a[8] * b[6];
+ r[7] = a[6] * b[1] + a[7] * b[4] + a[8] * b[7];
+ r[8] = a[6] * b[2] + a[7] * b[5] + a[8] * b[8];
+}
+
+/* computes conversion matrix from a given color space to CIE XYZ */
+void
+color_compute_color_space_to_xyz_matrix(double matrix[9], const struct color_space_chromacity_info *space)
+{
+ double wX = space->white[0] / space->white[1];
+ double wZ = (1 - space->white[0] - space->white[1]) / space->white[1];
+ double a[9], b[9];
+ a[0] = space->prim1[0]; a[3] = space->prim1[1]; a[6] = 1 - a[0] - a[3];
+ a[1] = space->prim2[0]; a[4] = space->prim2[1]; a[7] = 1 - a[1] - a[4];
+ a[2] = space->prim3[0]; a[5] = space->prim3[1]; a[8] = 1 - a[2] - a[5];
+ color_invert_matrix(b, a);
+ double ra = wX * b[0] + b[1] + wZ * b[2];
+ double rb = wX * b[3] + b[4] + wZ * b[5];
+ double rc = wX * b[6] + b[7] + wZ * b[8];
+ matrix[0] = a[0] * ra;
+ matrix[1] = a[1] * rb;
+ matrix[2] = a[2] * rc;
+ matrix[3] = a[3] * ra;
+ matrix[4] = a[4] * rb;
+ matrix[5] = a[5] * rc;
+ matrix[6] = a[6] * ra;
+ matrix[7] = a[7] * rb;
+ matrix[8] = a[8] * rc;
+}
+
+/* computes matrix to join transformations with different reference whites */
+void
+color_compute_bradford_matrix(double matrix[9], const double source[2], const double dest[2])
+{
+ /* cone response matrix and its inversion */
+ static double r[9] = {
+ 0.8951, 0.2664, -0.1614,
+ -0.7502, 1.7135, 0.0367,
+ 0.0389, -0.0685, 1.0296};
+ //static double i[9] = {0.9870, -0.1471, 0.1600, 0.4323, 0.5184, 0.0493, -0.0085, 0.0400, 0.9685};
+ double i[9];
+ color_invert_matrix(i, r);
+ double aX = source[0] / source[1];
+ double aZ = (1 - source[0] - source[1]) / source[1];
+ double bX = dest[0] / dest[1];
+ double bZ = (1 - dest[0] - dest[1]) / dest[1];
+ double x = (r[0] * bX + r[1] + r[2] * bZ) / (r[0] * aX + r[1] + r[2] * aZ);
+ double y = (r[3] * bX + r[4] + r[5] * bZ) / (r[3] * aX + r[4] + r[5] * aZ);
+ double z = (r[6] * bX + r[7] + r[8] * bZ) / (r[6] * aX + r[7] + r[8] * aZ);
+ double m[9];
+ m[0] = i[0] * x; m[1] = i[1] * y; m[2] = i[2] * z;
+ m[3] = i[3] * x; m[4] = i[4] * y; m[5] = i[5] * z;
+ m[6] = i[6] * x; m[7] = i[7] * y; m[8] = i[8] * z;
+ mul_matrices(matrix, m, r);
+}
+
+void
+color_compute_color_spaces_conversion_matrix(double matrix[9], const struct color_space_chromacity_info *src, const struct color_space_chromacity_info *dest)
+{
+ double a_to_xyz[9], b_to_xyz[9], xyz_to_b[9], bradford[9], m[9];
+ color_compute_color_space_to_xyz_matrix(a_to_xyz, src);
+ color_compute_color_space_to_xyz_matrix(b_to_xyz, dest);
+ color_invert_matrix(xyz_to_b, b_to_xyz);
+ if (src->white[0] == dest->white[0] && src->white[1] == dest->white[1])
+ mul_matrices(matrix, a_to_xyz, xyz_to_b);
+ else
+ {
+ color_compute_bradford_matrix(bradford, src->white, dest->white);
+ mul_matrices(m, a_to_xyz, bradford);
+ mul_matrices(matrix, m, xyz_to_b);
+ }
+}
+
+/* sRGB to XYZ */
+void
+srgb_to_xyz_exact(double xyz[3], double srgb[3])
+{
+ static double matrix[9] = {
+ 0.41248031, 0.35756952, 0.18043951,
+ 0.21268516, 0.71513904, 0.07217580,
+ 0.01933501, 0.11918984, 0.95031473};
+ double srgb_lin[3];
+ invert_gamma_detailed(srgb_lin, srgb, &color_srgb_info.gamma);
+ apply_matrix(xyz, srgb_lin, matrix);
+ xyz_to_srgb_exact(srgb_lin, xyz);
+}
+
+/* XYZ to sRGB */
+void
+xyz_to_srgb_exact(double srgb[3], double xyz[3])
+{
+ static double matrix[9] = {
+ 3.24026666, -1.53704957, -0.49850256,
+ -0.96928381, 1.87604525, 0.04155678,
+ 0.05564281, -0.20402363, 1.05721334};
+ double srgb_lin[3];
+ apply_matrix(srgb_lin, xyz, matrix);
+ clip(srgb_lin);
+ correct_gamma_detailed(srgb, srgb_lin, &color_srgb_info.gamma);