Add libred_get_colour_xy, libred_get_temperature_xy, and libred_get_temperature

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

6 files changed, 281 insertions, 43 deletions

diff --git a/Makefile b/Makefile
index 44116a7..ce4ce74 100644
--- a/Makefile
+++ b/Makefile
@@ -11,7 +11,7 @@ include mk/$(OS).mk
 
 
 LIB_MAJOR = 1
-LIB_MINOR = 0
+LIB_MINOR = 1
 LIB_VERSION = $(LIB_MAJOR).$(LIB_MINOR)
 
 
@@ -22,9 +22,13 @@ OBJ =\
 LOBJ = $(OBJ:.o=.lo)
 
 MAN0 = libred.h.0
-MAN3 = libred_get_colour.3 libred_solar_elevation.3
 MAN7 = libred.7
 
+MAN3 =\
+	libred_get_colour.3\
+	libred_get_colour_xy.3\
+	libred_solar_elevation.3
+
 
 all: libred.a libred.$(LIBEXT)
 solar.o: libred.h
@@ -66,9 +70,9 @@ install: libred.a libred.$(LIBEXT)
 	ln -sf -- libred.$(LIBMINOREXT) "$(DESTDIR)$(PREFIX)/lib/libred.$(LIBMAJOREXT)"
 	ln -sf -- libred.$(LIBMAJOREXT) "$(DESTDIR)$(PREFIX)/lib/libred.$(LIBEXT)"
 	cp -- libred.h "$(DESTDIR)$(PREFIX)/include"
-	cp -- $(MAN0) "$(DESTDIR)$(MANPREFIX)/man0"
-	cp -- $(MAN3) "$(DESTDIR)$(MANPREFIX)/man3"
-	cp -- $(MAN7) "$(DESTDIR)$(MANPREFIX)/man7"
+	cp -P -- $(MAN0) "$(DESTDIR)$(MANPREFIX)/man0"
+	cp -P -- $(MAN3) "$(DESTDIR)$(MANPREFIX)/man3"
+	cp -P -- $(MAN7) "$(DESTDIR)$(MANPREFIX)/man7"
 
 uninstall:
 	-rm -f -- "$(DESTDIR)$(PREFIX)/lib/libred.a"
diff --git a/blackbody.c b/blackbody.c
index 6edc85a..78608c1 100644
--- a/blackbody.c
+++ b/blackbody.c
@@ -81,7 +81,7 @@ ciexyy_to_srgb(double x, double y, double Y, double *r, double *g, double *b)
 
 /**
  * Perform linear interpolation (considered very good)
- * between the CIE xyY values for two colour temperatures
+ * between the CIE xy values for two colour temperatures
  * and convert the result to sRGB. The two colours should
  * be the closest below the desired colour temperature,
  * and the closest above the desired colour temperature
@@ -106,40 +106,29 @@ interpolate(double x1, double y1, double x2, double y2, double temp, double *r,
 }
 
 /**
- * Get the [0, 1] sRGB values of a colour temperature
- * 
- * libred has a table of colour temperature values, this
- * function interpolates values that are missing. If you
- * don't want any interpolation the `temp` parameter can
- * be specified in one of the following ways:
- * 
- * - floor:
- *         (temp - LIBRED_LOWEST_TEMPERATURE) /
- *         LIBRED_DELTA_TEMPERATURE *
- *         LIBRED_DELTA_TEMPERATURE +
- *         LIBRED_LOWEST_TEMPERATURE
- * 
- * - ceiling:
- *         (temp - LIBRED_LOWEST_TEMPERATURE +
- *          LIBRED_DELTA_TEMPERATURE - 1) /
- *         LIBRED_DELTA_TEMPERATURE *
- *         LIBRED_DELTA_TEMPERATURE +
- *         LIBRED_LOWEST_TEMPERATURE
- * 
- * - round to nearest:
- *         (temp - LIBRED_LOWEST_TEMPERATURE +
- *          LIBRED_DELTA_TEMPERATURE / 2) /
- *         LIBRED_DELTA_TEMPERATURE *
- *         LIBRED_DELTA_TEMPERATURE +
- *         LIBRED_LOWEST_TEMPERATURE
+ * Perform linear interpolation (considered very good)
+ * between the CIE xy values for two colour temperatures.
+ * The two colours should be the closest below the desired
+ * colour temperature, and the closest above the desired
+ * colour temperature
  * 
- * @param   temp  The desired colour temperature
- * @param   r     Output parameter for the “red” value
- * @param   g     Output parameter for the green value
- * @param   b     Output parameter for the blue value
- * @return        0 on succeess, -1 on error
- * @throws  EDOM  The selected temperature is below 1000K
+ * @param  x1    The 'x' component for the low colour
+ * @param  y1    The 'y' component for the low colour
+ * @param  x2    The 'x' component for the high colour
+ * @param  y2    The 'y' component for the high colour
+ * @param  temp  The desired colour temperature
+ * @param  x     Output parameter for the CIE x value
+ * @param  y     Output parameter for the CIE y value
  */
+static void
+interpolate_xy(double x1, double y1, double x2, double y2, double temp, double *x, double *y)
+{
+	double weight = fmod(temp - (LIBRED_LOWEST_TEMPERATURE % LIBRED_DELTA_TEMPERATURE),
+	                     (double)LIBRED_DELTA_TEMPERATURE) / (double)LIBRED_DELTA_TEMPERATURE;
+	*x = x1 * (1 - weight) + x2 * weight;
+	*y = y1 * (1 - weight) + y2 * weight;
+}
+
 int
 libred_get_colour(long int temp, double *r, double *g, double *b)
 {
@@ -156,7 +145,7 @@ libred_get_colour(long int temp, double *r, double *g, double *b)
 	}
 
 	tmp = temp - LIBRED_LOWEST_TEMPERATURE;
-	
+
 	i = (size_t)(tmp / LIBRED_DELTA_TEMPERATURE);
 	if (tmp % LIBRED_DELTA_TEMPERATURE) {
 		x1 = xy_table[i].x;
@@ -172,6 +161,139 @@ libred_get_colour(long int temp, double *r, double *g, double *b)
 
 	return 0;
 }
+
+int
+libred_get_colour_xy(long int temp, double *x, double *y) /* TODO test */
+{
+	double x1, y1, x2, y2;
+	size_t i;
+	long int tmp;
+
+	if (temp > LIBRED_HIGHEST_TEMPERATURE)
+		temp = LIBRED_HIGHEST_TEMPERATURE;
+
+	if (temp < LIBRED_LOWEST_TEMPERATURE) {
+		errno = EDOM;
+		return -1;
+	}
+
+	tmp = temp - LIBRED_LOWEST_TEMPERATURE;
+
+	i = (size_t)(tmp / LIBRED_DELTA_TEMPERATURE);
+	if (tmp % LIBRED_DELTA_TEMPERATURE) {
+		x1 = xy_table[i].x;
+		y1 = xy_table[i].y;
+		x2 = xy_table[i + 1].x;
+		y2 = xy_table[i + 1].y;
+		interpolate_xy(x1, y1, x2, y2, (double)temp, x, y);
+	} else {
+		*x = xy_table[i].x;
+		*y = xy_table[i].y;
+	}
+
+	return 0;
+}
+
+double
+libred_get_temperature_xy(double x, double y, double *x_error, double *y_error) /* TODO man, test */
+{
+	size_t i, j;
+	double x1, y1, x2, y2, dx, dy, xd, yd, t, d2;
+	double best_temp = 0, best_d2 = INFINITY;
+
+	if (!x_error)
+		x_error = &dx;
+	if (!y_error)
+		y_error = &dy;
+
+	x1 = xy_table[0].x;
+	y1 = xy_table[0].y;
+
+	for (j = 1; j < sizeof(xy_table) / sizeof(*xy_table); j++) {
+		x2 = xy_table[j].x;
+		y2 = xy_table[j].y;
+		dx = x2 - x1;
+		dy = y2 - y1;
+
+		t = ((x - x1) * dx + (y - y1) * dy) / (dx * dx + dy * dy);
+		if (!isfinite(t))
+			continue;
+		t = t < 0 ? 0 : t;
+		t = t > 1 ? 1 : t;
+
+		xd = dx * t + x1 - x;
+		yd = dy * t + y1 - y;
+		d2 = xd * xd + yd * yd;
+
+		if (d2 < best_d2) {
+			*x_error = xd;
+			*y_error = yd;
+			best_d2 = d2;
+			t *= (double)(j - i);
+			best_temp = (double)i + t;
+		}
+	}
+
+	return best_temp * LIBRED_DELTA_TEMPERATURE + LIBRED_LOWEST_TEMPERATURE;
+}
+
+/**
+ * Convert from [0, 1] sRGB to CIE xyY
+ * 
+ * @param  r  The “red” value
+ *            (Seriously, sRGB red is orange, just look at it fullscreen)
+ * @param  g  The green value
+ * @param  b  The blue value
+ * @param  x  Output parameter for The 'x' component
+ * @param  y  Output parameter for The 'y' component
+ * @param  Y  Output parameter for The 'Y' component
+ */
+static void
+srgb_to_ciexyy(double r, double g, double b, double *x, double *y, double *Y)
+{
+	double s, z;
+
+	r *= s = r < 0 ? -1 : 1;
+	r = s * (r <= 0.0031306684425217108 * 12.92 ? r * 12.92 : pow((r + 0.055) / 1.055, 2.4));
+	g *= s = g < 0 ? -1 : 1;
+	g = s * (g <= 0.0031306684425217108 * 12.92 ? g * 12.92 : pow((g + 0.055) / 1.055, 2.4));
+	b *= s = b < 0 ? -1 : 1;
+	b = s * (g <= 0.0031306684425217108 * 12.92 ? b * 12.92 : pow((b + 0.055) / 1.055, 2.4));
+
+	*x = (0.41245744558236758 * r) + (0.35757586524551588 * g) + (0.18043724782639967 * b);
+	*y = (0.21267337037840828 * r) + (0.71515173049103176 * g) + (0.07217489913055987 * b);
+	z  = (0.01933394276167346 * r) + (0.11919195508183859 * g) + (0.95030283855237174 * b);
+
+	*Y = *y;
+	s = *x + *y + z;
+	*x /= s;
+	*y /= s;
+	if (!isfinite(*x) || !isfinite(*y))
+		*x = *y = 0;
+}
+
+double
+libred_get_temperature(double r, double g, double b, double *y, /* TODO man, test */
+                       double *r_error, double *g_error, double *b_error)
+{
+	double tx, ty, luma, x_error, y_error, ret;
+
+	srgb_to_ciexyy(r, g, b, &tx, &ty, &luma);
+	ret = libred_get_temperature_xy(tx, ty, &x_error, &y_error);
+	ciexyy_to_srgb(x_error, y_error, luma, &r, &g, &b);
+
+	if (y)
+		*y = luma;
+	if (r_error)
+		*r_error = r;
+	if (g_error)
+		*g_error = g;
+	if (b_error)
+		*b_error = b;
+
+	return ret;
+}
+
 #endif
 
 
diff --git a/libred.h b/libred.h
index 3b0e489..29a00cd 100644
--- a/libred.h
+++ b/libred.h
@@ -212,6 +212,82 @@ int libred_check_timetravel(void);
  */
 int libred_get_colour(long int, double *, double *, double *);
 
+/**
+ * Get the CIE xy values of a colour temperature
+ * 
+ * libred has a table of colour temperature values, this
+ * function interpolates values that are missing. If you
+ * don't want any interpolation the `temp` parameter can
+ * be specified in one of the following ways:
+ * 
+ * - floor:
+ *         (temp - LIBRED_LOWEST_TEMPERATURE) /
+ *         LIBRED_DELTA_TEMPERATURE *
+ *         LIBRED_DELTA_TEMPERATURE +
+ *         LIBRED_LOWEST_TEMPERATURE
+ * 
+ * - ceiling:
+ *         (temp - LIBRED_LOWEST_TEMPERATURE +
+ *          LIBRED_DELTA_TEMPERATURE - 1) /
+ *         LIBRED_DELTA_TEMPERATURE *
+ *         LIBRED_DELTA_TEMPERATURE +
+ *         LIBRED_LOWEST_TEMPERATURE
+ * 
+ * - round to nearest:
+ *         (temp - LIBRED_LOWEST_TEMPERATURE +
+ *          LIBRED_DELTA_TEMPERATURE / 2) /
+ *         LIBRED_DELTA_TEMPERATURE *
+ *         LIBRED_DELTA_TEMPERATURE +
+ *         LIBRED_LOWEST_TEMPERATURE
+ * 
+ * @param   temp  The desired colour temperature
+ * @param   x     Output parameter for the x value
+ * @param   y     Output parameter for the y value
+ * @return        0 on succeess, -1 on error
+ * 
+ * @throws  EDOM  The selected temperature is below 1000K
+ */
+int libred_get_colour_xy(long int, double *, double *);
+
+/**
+ * Calculate the colour temperature from its CIE xy values
+ * 
+ * @param   x        The CIE x value for the colour temperature
+ * @param   y        The CIE y value for the colour temperature
+ * @param   x_error  Output parameter for the absolute error in
+ *                   the CIE x value of the returned colour
+ *                   temperature; may be `NULL`
+ * @param   y_error  Output parameter for the absolute error in
+ *                   the CIE y value of the returned colour
+ *                   temperature; may be `NULL`
+ * @return           The closest matching colour temperature
+ */
+double libred_get_temperature_xy(double, double, double *, double *);
+
+/**
+ * Calculate the colour temperature from its [0, 1] sRGB values
+ * 
+ * @param   r        The sRGB “red” value for the colour temperature
+ * @param   g        The sRGB green value for the colour temperature
+ * @param   b        The sRGB blue value for the colour temperature
+ * @param   y        Output parameter for the CIE Y value for the
+ *                   provided sRGB colour, the sRGB error values must
+ *                   be multiplied by this value to correct them to
+ *                   match the brightness of the input sRGB colour;
+ *                   may be `NULL`
+ * @param   r_error  Output parameter for the absolute error in the
+ *                   “red” value of the returned colour temperature;
+ *                   may be `NULL`
+ * @param   g_error  Output parameter for the absolute error in the
+ *                   green value of the returned colour temperature;
+ *                   may be `NULL`
+ * @param   b_error  Output parameter for the absolute error in the
+ *                   blue value of the returned colour temperature;
+ *                   may be `NULL`
+ * @return           The closest matching colour temperature
+ */
+double libred_get_temperature(double, double, double, double *, double *, double *, double *);
+
 
 
 #endif
diff --git a/libred_get_colour.3 b/libred_get_colour.3
index afd0a7f..61a62ff 100644
--- a/libred_get_colour.3
+++ b/libred_get_colour.3
@@ -6,6 +6,7 @@ libred_get_colour \- Calculate a colour temperature
 #include <libred.h>
 
 int \fBlibred_get_colour\fP(long int \fItemp\fP, double *\fIr\fP, double *\fIg\fP, double *\fIb\fP);
+int \fBlibred_get_colour_xy\fP(long int \fItemp\fP, double *\fIx\fP, double *\fIy\fP);
 .fi
 .PP
 Link with
@@ -14,7 +15,8 @@ Link with
 .BR libred_get_colour ()
 gets or interpolates the colour temperature for
 .I temp
-kelvins, and returns the colour temperature in sRGB. The values,
+kelvins, and returns the colour temperature in sRGB, with the
+sRGB transfer function applied. The values,
 between 0.0 and 1.0, for the \(dqred\(dq, green, and blue channels
 are stored in
 .IR *r ,
@@ -49,14 +51,45 @@ ceiling
 .TP
 round to nearest
 .I (temp - LIBRED_LOWEST_TEMPERATURE + LIBRED_DELTA_TEMPERATURE / 2) / LIBRED_DELTA_TEMPERATURE * LIBRED_DELTA_TEMPERATURE + LIBRED_LOWEST_TEMPERATURE
+.PP
+The
+.BR libred_get_colour_xy ()
+function is the same as the
+.BR libred_get_colour ()
+except it uses the CIE xy colour space instead of sRGB
+(of course without any transfer function). This means
+that the chroma values x and y are stored in
+.I *x
+and
+.IR *y ,
+and the CIE Y value shall be assumed to be 1 for
+equivalence with the
+.BR libred_get_colour ()
+function.
+.PP
+Neither
+.IR r ,
+.IR g ,
+.IR b ,
+.IR x ,
+nor
+.I y
+may be
+.IR NULL .
 .SH RETURN VALUE
 Upon successful completion, the
 .BR libred_get_colour ()
+and
+.BR libred_get_colour_xy ()
 function returns 0. On failure, the function returns -1 and sets
 .I errno
 to indicate the error.
 .SH ERRORS
-The function may fail if:
+The
+.BR libred_get_colour ()
+and
+.BR libred_get_colour_xy ()
+functions may fail if:
 .TP
 .B EDOM
 If
diff --git a/libred_get_colour_xy.3 b/libred_get_colour_xy.3
new file mode 120000
index 0000000..1d7ee17
--- /dev/null
+++ b/libred_get_colour_xy.3
@@ -0,0 +1 @@
+libred_get_colour.3
\ No newline at end of file
diff --git a/libred_solar_elevation.3 b/libred_solar_elevation.3
index 9bb4b0a..70587ca 100644
--- a/libred_solar_elevation.3
+++ b/libred_solar_elevation.3
@@ -27,11 +27,13 @@ Other values may or may not work, no error is thrown if used.
 .SH RETURN VALUE
 Upon successful completion, the 
 .BR libred_solar_elevation ()
-return 0, on failure it returns -1 and sets
+function returns 0, on failure it returns -1 and sets
 .I errno
 to indicate the error.
 .SH ERRORS
-The function may fail for any reason specified for
+The
+.BR libred_solar_elevation ()
+function may fail for any reason specified for
 .BR clock_gettime (3).
 .SH SEE ALSO
 .BR libred.h (0),