Add RGB colourspace conversion

Signed-off-by: Mattias Andrée <maandree@kth.se>

2 files changed, 359 insertions, 0 deletions

diff --git a/src/libclut.c b/src/libclut.c
index 8963a4d..93d0ac0 100644
--- a/src/libclut.c
+++ b/src/libclut.c
@@ -16,6 +16,8 @@
  */
 #include "libclut.h"
 
+#include <errno.h>
+
 
 
 /**
@@ -199,3 +201,215 @@ void (libclut_model_cielab_to_ciexyz)(double L, double a, double b, double* X, d
   libclut_model_cielab_to_ciexyz(L, a, b, X, Y, Z);
 }
 
+
+
+/**
+ * Divide all values in a row by a divisor.
+ * 
+ * @param  m  The first part of the row.
+ * @param  a  The second part of the row.
+ * @param  d  The divisor.
+ */
+static void divrow(double m[3], double a[3], double d)
+{
+  m[0] /= d, m[1] /= d, m[2] /= d;
+  a[0] /= d, a[1] /= d, a[2] /= d;
+}
+
+/**
+ * Subtract all values in a row by corresponding value
+ * in another row multiplied by a common value.
+ * 
+ * @param  a1  The first part of the minuend/difference row.
+ * @param  a2  The second part of the minuend/difference row.
+ * @param  b1  The first part of the subtrahend row.
+ * @param  b2  The second part of the subtrahend row.
+ * @param  m   The multiplier.
+ */
+static void subrow(double a1[3], double a2[3], double b1[3], double b2[3], double m)
+{
+  a1[0] -= b1[0] * m, a1[1] -= b1[1] * m, a1[2] -= b1[2] * m;
+  a2[0] -= b2[0] * m, a2[1] -= b2[1] * m, a2[2] -= b2[2] * m;
+}
+
+
+/**
+ * Invert a matrix.
+ * 
+ * @param   M  The matrix to invert, will be modified to an
+ *             identity matrix, possibly with reordered rows.
+ * @param   A  The inversion of M (as input).
+ * @return     1 on success, 0 if the matrix is not invertible.
+ */
+static int invert(libclut_colourspace_conversion_matrix_t M, libclut_colourspace_conversion_matrix_t A)
+{
+  int r0 = 0, r1 = 1, r2 = 2, t, swapped = 0;
+  libclut_colourspace_conversion_matrix_t T;
+  
+  A[0][0] = A[1][1] = A[2][2] = 1;
+  A[0][1] = A[0][2] = A[1][0] = A[1][2] = A[2][0] = A[2][1] = 0;
+  
+  if (libclut_0__(M[r0][0]))
+    {
+      if (libclut_0__(M[r1][0]))
+	{
+	  if (libclut_0__(M[r2][0]))
+	    return 0;
+	  t = r0, r0 = r2, r2 = t;
+	}
+      else
+	t = r0, r0 = r1, r1 = t;
+      swapped = 1;
+    }
+  
+  divrow(M[r0], A[r0], M[r0][0]);
+  subrow(M[r1], A[r1], M[r0], A[r0], M[r1][0]);
+  subrow(M[r2], A[r2], M[r0], A[r0], M[r2][0]);
+  
+  if (libclut_0__(M[r1][1]))
+    {
+      if (libclut_0__(M[r2][1]))
+	return 0;
+      t = r1, r1 = r2, r2 = t;
+      swapped = 1;
+    }
+  
+  divrow(M[r1], A[r1], M[r1][1]);
+  subrow(M[r2], A[r2], M[r1], A[r1], M[r2][1]);
+  
+  if (libclut_0__(M[r2][2]))
+    return 0;
+  
+  divrow(M[r2], A[r2], M[r2][2]);
+  
+  subrow(M[r1], A[r1], M[r2], A[r2], M[r1][2]);
+  subrow(M[r0], A[r0], M[r2], A[r2], M[r0][2]);
+  
+  subrow(M[r0], A[r0], M[r1], A[r1], M[r0][1]);
+  
+  if (swapped)
+    {
+      memcpy(T, A, sizeof(T));
+      memcpy(A[0], T[r0], sizeof(*T));
+      memcpy(A[1], T[r1], sizeof(*T));
+      memcpy(A[2], T[r2], sizeof(*T));
+    }
+  
+  return 1;
+}
+
+
+/**
+ * Create an RGB to CIE XYZ conversion matrix.
+ * 
+ * @param   cs  The colour space.
+ * @param   M   The output matrix.
+ * @return      Zero on success, -1 on error.
+ * 
+ * @throws  EINVAL  The colourspace cannot be used.
+ */
+static int get_conversion_matrix(const libclut_rgb_colourspace_t* cs, libclut_colourspace_conversion_matrix_t M)
+{
+#define XYY_TO_XYZ(x, y, Y, Xp, Yp, Zp)		\
+  (libclut_0__(Y)) ?				\
+    (*(Xp) = *(Zp) = 0, *(Yp) = (Y)) :		\
+    (*(Xp) = (x) * (Y) / (y),			\
+     *(Yp) = (Y),				\
+     *(Zp) = (1 - (x) - (y)) * (Y) / (y))
+  
+  double Xr, Yr, Zr, Xg, Yg, Zg, Xb, Yb, Zb, Xw, Yw, Zw, Sr, Sg, Sb;
+  libclut_colourspace_conversion_matrix_t M2;
+  
+  XYY_TO_XYZ(cs->red_x,   cs->red_y,   1,           &Xr, &Yr, &Zr);
+  XYY_TO_XYZ(cs->green_x, cs->green_y, 1,           &Xg, &Yg, &Zg);
+  XYY_TO_XYZ(cs->blue_x,  cs->blue_y,  1,           &Xb, &Yb, &Zb);
+  XYY_TO_XYZ(cs->white_x, cs->white_y, cs->white_Y, &Xw, &Yw, &Zw);
+  
+  M2[0][0] = Xr, M2[0][1] = Xg, M2[0][2] = Xb;
+  M2[1][0] = Yr, M2[1][1] = Yg, M2[1][2] = Yb;
+  M2[2][0] = Zr, M2[2][1] = Zg, M2[2][2] = Zb;
+  
+  if (!invert(M2, M))
+    return errno = EINVAL, -1;
+  
+  Sr = M[0][0] * Xw + M[0][1] * Yw + M[0][2] * Zw;
+  Sg = M[1][0] * Xw + M[1][1] * Yw + M[1][2] * Zw;
+  Sb = M[2][0] * Xw + M[2][1] * Yw + M[2][2] * Zw;
+  
+  M[0][0] = Sr * Xr, M[0][1] = Sg * Xg, M[0][2] = Sb * Xb;
+  M[1][0] = Sr * Yr, M[1][1] = Sg * Yg, M[1][2] = Sb * Yb;
+  M[2][0] = Sr * Zr, M[2][1] = Sg * Zg, M[2][2] = Sb * Zb;
+  
+  return 0;
+  
+#undef XYY_TO_XYZ
+}
+
+
+/**
+ * Create a matrix for converting values between
+ * two RGB colourspaces.
+ * 
+ * @param   from  The input colourspace, the Y-component is only necessary for the whitepoint.
+ * @param   to    The output colourspace, the Y-component is only necessary for the whitepoint.
+ * @param   M     Output matrix for conversion from `from` to `to`.
+ * @param   Minv  Output matrix for conversion from `to` to `from`, may be `NULL`.
+ * @return        Zero on success, -1 on error.
+ * 
+ * @throws  EINVAL  The colourspace cannot be used.
+ */
+int libclut_model_get_rgb_conversion_matrix(const libclut_rgb_colourspace_t* from,
+					    const libclut_rgb_colourspace_t* to,
+					    libclut_colourspace_conversion_matrix_t M,
+					    libclut_colourspace_conversion_matrix_t Minv)
+{
+  libclut_colourspace_conversion_matrix_t A, B;
+
+  if (get_conversion_matrix(from, A))
+    return -1;
+  if (get_conversion_matrix(to, M))
+    return -1;
+  if (!invert(M, B))
+    return errno = EINVAL, -1;
+  
+  M[0][0] = A[0][0] * B[0][0] + A[0][1] * B[1][0] + A[0][2] * B[2][0];
+  M[0][1] = A[0][0] * B[0][1] + A[0][1] * B[1][1] + A[0][2] * B[2][1];
+  M[0][2] = A[0][0] * B[0][2] + A[0][1] * B[1][2] + A[0][2] * B[2][2];
+  
+  M[1][0] = A[1][0] * B[0][0] + A[1][1] * B[1][0] + A[1][2] * B[2][0];
+  M[1][1] = A[1][0] * B[0][1] + A[1][1] * B[1][1] + A[1][2] * B[2][1];
+  M[1][2] = A[1][0] * B[0][2] + A[1][1] * B[1][2] + A[1][2] * B[2][2];
+  
+  M[2][0] = A[2][0] * B[0][0] + A[2][1] * B[1][0] + A[2][2] * B[2][0];
+  M[2][1] = A[2][0] * B[0][1] + A[2][1] * B[1][1] + A[2][2] * B[2][1];
+  M[2][2] = A[2][0] * B[0][2] + A[2][1] * B[1][2] + A[2][2] * B[2][2];
+  
+  if (Minv != NULL)
+    {
+      memcpy(A, M, sizeof(A));
+      if (!invert(A, Minv))
+	return errno = EINVAL, -1;
+    }
+  
+  return 0;
+}
+
+
+/**
+ * Convert an RGB colour into another RGB colourspace.
+ * None of the parameter may have side-effects.
+ * 
+ * @param  r      The red component of the colour to convert.
+ * @param  g      The green component of the colour to convert.
+ * @param  b      The blue component of the colour to convert.
+ * @param  M      Conversion matrix, create with `libclut_model_get_rgb_conversion_matrix`.
+ * @param  out_r  Output parameter for the new red component.
+ * @param  out_g  Output parameter for the new green component.
+ * @param  out_b  Output parameter for the new blue component.
+ */
+void (libclut_model_convert_rgb)(double r, double g, double b, libclut_colourspace_conversion_matrix_t M,
+				  double *out_r, double *out_g, double *out_b)
+{
+  libclut_model_convert_rgb(r, g, b, M, out_r, out_g, out_b);
+}
+
diff --git a/src/libclut.h b/src/libclut.h
index 973ac0a..9f3326b 100644
--- a/src/libclut.h
+++ b/src/libclut.h
@@ -23,6 +23,94 @@
 
 
 
+/**
+ * Initialiser for `struct libclut_rgb_colourspace` with the values
+ * of the sRGB colour space.
+ */
+#define LIBCLUT_RGB_COLOURSPACE_SRGB_INITIALISER			\
+  {									\
+    .red_x   = 0.6400,   .red_y   = 0.3300,   .red_Y   = 0.212656,	\
+    .green_x = 0.3000,   .green_y = 0.6000,   .green_Y = 0.715158,	\
+    .blue_x  = 0.1500,   .blue_y  = 0.0600,   .blue_Y  = 0.072186,	\
+    .white_x = 0.31271,  .white_y = 0.32902,  .white_Y = 1.0000		\
+  }
+
+
+
+/**
+ * RGB colour space structure.
+ */
+typedef struct libclut_rgb_colourspace
+{
+  /**
+   * The x-component of the red colour's xyY value.
+   */
+  double red_x;
+  
+  /**
+   * The y-component of the red colour's xyY value.
+   */
+  double red_y;
+  
+  /**
+   * The Y-component of the red colour's xyY value.
+   */
+  double red_Y;
+  
+  /**
+   * The x-component of the green colour's xyY value.
+   */
+  double green_x;
+  
+  /**
+   * The y-component of the green colour's xyY value.
+   */
+  double green_y;
+  
+  /**
+   * The Y-component of the green colour's xyY value.
+   */
+  double green_Y;
+  
+  /**
+   * The x-component of the blue colour's xyY value.
+   */
+  double blue_x;
+  
+  /**
+   * The y-component of the blue colour's xyY value.
+   */
+  double blue_y;
+  
+  /**
+   * The Y-component of the blue colour's xyY value.
+   */
+  double blue_Y;
+  
+  /**
+   * The x-component of the white point's xyY value.
+   */
+  double white_x;
+  
+  /**
+   * The y-component of the white point's xyY value.
+   */
+  double white_y;
+  
+  /**
+   * The Y-component of the white point's xyY value.
+   */
+  double white_Y;
+} libclut_rgb_colourspace_t;
+
+
+/**
+ * Matrix date-type for colourspace conversion.
+ */
+typedef double libclut_colourspace_conversion_matrix_t[3][3];
+
+
+
 /* This is to avoid warnings about comparing double, These are only
  * used when it is safe, for example to test whether optimisations
  * are possible. { */
@@ -1082,6 +1170,9 @@ static inline int libclut_0__(double x)  {  return libclut_eq__(x, 0);  }
   (size_t)((double)(i) * (double)(out) / (double)(in))
 
 
+/* TODO libclut_convert_rgb */
+
+
 
 #if defined(__GNUC__) && !defined(__clang__)
 # define LIBCLUT_GCC_ONLY__(x)  x
@@ -1390,6 +1481,60 @@ void (libclut_model_cielab_to_ciexyz)(double, double, double, double*, double*,
   (((C)*(C)*(C) > 0.00885642) ? ((C)*(C)*(C)) : (((C) - 0.1379310) / (7.78 + 703.0 / 99900)))
 
 
+/**
+ * Create a matrix for converting values between
+ * two RGB colourspaces.
+ * 
+ * @param   from  The input colourspace, the Y-component is only necessary for the whitepoint.
+ * @param   to    The output colourspace, the Y-component is only necessary for the whitepoint.
+ * @param   M     Output matrix for conversion from `from` to `to`.
+ * @param   Minv  Output matrix for conversion from `to` to `from`, may be `NULL`.
+ * @return        Zero on success, -1 on error.
+ * 
+ * @throws  EINVAL  The colourspace cannot be used.
+ */
+LIBCLUT_GCC_ONLY__(__attribute__((__leaf__)))
+int libclut_model_get_rgb_conversion_matrix(const libclut_rgb_colourspace_t*,
+					    const libclut_rgb_colourspace_t*,
+					    libclut_colourspace_conversion_matrix_t,
+					    libclut_colourspace_conversion_matrix_t);
+/* TODO doc libclut_model_get_rgb_conversion_matrix */
+
+
+/**
+ * Convert an RGB colour into another RGB colourspace.
+ * None of the parameter may have side-effects.
+ * 
+ * Requires linking with '-lclut', or '-lm' if
+ * `libclut_model_standard_to_linear1` or
+ * `libclut_model_linear_to_standard1` is not undefined.
+ * 
+ * @param  r      The red component of the colour to convert.
+ * @param  g      The green component of the colour to convert.
+ * @param  b      The blue component of the colour to convert.
+ * @param  M      Conversion matrix, create with `libclut_model_get_rgb_conversion_matrix`,
+ *                must not have side-effects (unless libclut_model_convert_rgb1 is undefined).
+ * @param  out_r  Output parameter for the new red component.
+ * @param  out_g  Output parameter for the new green component.
+ * @param  out_b  Output parameter for the new blue component.
+ */
+LIBCLUT_GCC_ONLY__(__attribute__((__leaf__)))
+void (libclut_model_convert_rgb)(double, double, double, libclut_colourspace_conversion_matrix_t,
+				  double *, double *, double *);
+#define libclut_model_convert_rgb(r, g, b, M, out_r, out_g, out_b)						\
+  do														\
+    {														\
+      double r__ = libclut_model_standard_to_linear1(r);							\
+      double g__ = libclut_model_standard_to_linear1(g);							\
+      double b__ = libclut_model_standard_to_linear1(b);							\
+      *(out_r) = libclut_model_linear_to_standard1((M)[0][0] * r__ + (M)[0][1] * g__ + (M)[0][2] * b__);	\
+      *(out_g) = libclut_model_linear_to_standard1((M)[1][0] * r__ + (M)[1][1] * g__ + (M)[1][2] * b__);	\
+      *(out_b) = libclut_model_linear_to_standard1((M)[2][0] * r__ + (M)[2][1] * g__ + (M)[2][2] * b__);	\
+    }														\
+  while (0)
+/* TODO libclut_model_convert_rgb */
+
+
 
 #if defined(__clang__)
 # pragma GCC diagnostic pop