summaryrefslogtreecommitdiffstats
path: root/color.c
blob: ce34b5bb763ecc3feb45bc220f5639ee1085e8b7 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
/*
 * Color space conversion functions, as well as functions to compute
 * color coordinates of locations of the spectral and Planck (thermal)
 * loci.  The color space used is CIE xy; combine with intensity (Y)
 * for a total set of color coordinates.  Colors need to be converted
 * to XYZ space (weighted by intensity) before mixing.
 */

#include "color.h"

/* Convert CIE XYZ to CIE xy */
struct cie_xy const_func XYZ_to_xy(struct cie_XYZ XYZ)
{
	struct cie_xy xy;
	double total = XYZ.X + XYZ.Y + XYZ.Z;

	xy.x = XYZ.X / total;
	xy.y = XYZ.Y / total;

	return xy;
}

/* Convert CIE xy + Y to CIE XYZ */
struct cie_XYZ const_func xyY_to_XYZ(struct cie_xy xy, double Y)
{
	struct cie_XYZ XYZ;
	double Y_y = Y / xy.y;	/* Intensity weighting factor */

	XYZ.Y = Y;
	XYZ.X = Y_y * xy.x;
	XYZ.Z = Y_y - Y_y * (xy.x + xy.y);

	return XYZ;
}

struct Luv_uv_prime {
	double u, v;		/* Intermediate u' and v' values */
};

static struct const_func Luv_uv_prime Luv_uv_prime(struct cie_XYZ XYZ)
{
	struct Luv_uv_prime uv;
	double den = XYZ.X + 15.0*XYZ.Y + 3.0*XYZ.Z;

	uv.u = 4.0*XYZ.X / den;
	uv.v = 9.0*XYZ.Y / den;

	return uv;
}

/* Convert CIE XYZ to CIE L*u*v* given a specific whitepoint */
struct cie_Luv const_func XYZ_to_Luv(struct cie_XYZ XYZ, struct cie_XYZ white)
{
	struct cie_Luv Luv;
	struct Luv_uv_prime uv;	       /* u'v' */
	struct Luv_uv_prime uv_n;      /* u'v' for whitepoint */
	double Y_Yn = XYZ.Y / white.Y; /* Relative intensity */
	const double Y_cutoff = (6.0*6.0*6.0)/(29.0*29.0*29.0);
	double L13;

	if (Y_Yn <= Y_cutoff)
		Luv.L = ((29.0*29.0*29.0)/(3.0*3.0*3.0)) * Y_Yn;
	else
		Luv.L = 116.0*cbrt(Y_Yn) - 16.0;

	L13 = 13.0 * Luv.L;

	uv   = Luv_uv_prime(XYZ);
	uv_n = Luv_uv_prime(white);

	Luv.u = L13 * (uv.u - uv_n.u);
	Luv.v = L13 * (uv.v - uv_n.v);

	return Luv;
}

/* Convert CIE L*u*v* to CIE XYZ */
struct cie_XYZ const_func Luv_to_XYZ(struct cie_Luv Luv, struct cie_XYZ white)
{
	struct Luv_uv_prime uv_n; /* u'v' for white point */
	double u, v;		  /* u'v' */
	struct cie_XYZ XYZ;
	double Y;
	double L13_inv;		/* 1.0/(13.0*L) */

	uv_n = Luv_uv_prime(white);

	L13_inv = 1.0/(13.0 * Luv.L);

	u = Luv.u * L13_inv + uv_n.u;
	v = Luv.v * L13_inv + uv_n.v;

	if (Luv.L <= 8.0) {
		Y = white.Y * (Luv.L * ((3.0*3.0*3.0)/(29.0*29.0*29.0)));
	} else {
		double q = (Luv.L + 16.0)/116.0;
		Y = white.Y * (q * q * q);
	}

	XYZ.X = 9.0 * Y * u / (4.0 * v);
	XYZ.Y = Y;
	XYZ.Z = - Y * (3.0*u + 20.0*v - 12.0)/(4.0*v);

	return XYZ;
}

/*
 * Compute CIE xy for a specific spectral wavelength in nm.
 * The Y coordinates can be obtained from the photometric flux
 * (given in candela.)
 *
 * These are the 10° CIE 2006 values.
 */
#include "lin2012xyz10e_fine_7sf.c"

struct cie_xy pure_func spectral(double nm)
{
	int tblstep = round(nm * xytbl_step);

	if (tblstep < xytbl_first)
		return xytbl[0]; /* Arbitrary */
	else if (tblstep > xytbl_last)
		return xytbl[xytbl_last - xytbl_first]; /* Arbitrary */
	else
		return xytbl[tblstep - xytbl_first];
}

/*
 * Return the photometric scaling factor for monochromatic light of a
 * specific wavelength.  This factor applies:
 *
 * <photometric> = <bolometric>  * <factor>
 * <bolometric>  = <photometric> / <factor>
 *
 * This applies to the following units/quantities:
 *
 * Quantity		Bolometric unit			Photometric unit
 * ---------------------------------------------------------------------
 * Intensity		W/sr (watt per steradian)	cd (candela)
 * Luminous flux	W (watt)			lm (lumen)
 * Illuminance		W/m² (watt per square meter)	lx (lux)
 */
double photometric_scale(double nm)
{
	int tblstep = round(nm * xytbl_step);

	if (tblstep < xytbl_first || tblstep > xytbl_last)
		return 0.0;
	else
		return XYZtbl[tblstep - xytbl_first].Y * 683.0;
}

/*
 * Compute CIE XYZ coordinates for a particular point on the Planckian
 * locus (color temperature)
 */
struct cie_XYZ pure_func planckian(double T)
{
	/*
	 * Second radiation constant c2 = hc/k, value from CODATA:
	 * https://physics.nist.gov/cgi-bin/cuu/Value?c22ndrc|search_for=radiation
	 */
	const double c2 = 1.43877736e-2; /* m·K */
	double M;
	int i;
	struct cie_XYZ XYZ;
	const struct cie_XYZ *p;

	XYZ.X = XYZ.Y = XYZ.Z = 0.0;

	p = &XYZtbl[0];
	for (i = xytbl_first; i <= xytbl_last; i++) {
		double lambda = i * xytbl_m;
		double lambda2, lambda5;

		lambda2 = lambda*lambda;
		lambda5 = lambda2*lambda2*lambda;

		M = 1.0/(lambda5 * expm1(c2/(lambda*T)));
		XYZ.X += p->X * M;
		XYZ.Y += p->Y * M;
		XYZ.Z += p->Z * M;

		p++;
	}

	return XYZ;
}

/*
 * Convert sRGB color space, the default color space for image files,
 * to/from XYZ.  This is VERY different color space from the LED RGB
 * color spaces: not only does it have a much smaller color gamut, but
 * it is also nonlinear.
 */
static double const_func sRGB_to_linear(double c)
{
	if (c <= 0.04045)
		c = c/(255.0*12.92);
	else
		c = pow((c + 0.055*255.0)/(255.0*1.055), 2.4);

	return c;
}

static double const_func linear_to_sRGB(double c)
{
	/* Cap the value to the range [0,1] */
	c = signbit(c) ? 0.0 : (c > 1.0) ? 1.0 : c;

	if (c <= 1.0/0.04045)
		c = (12.92*255.0) * c;
	else
		c = (255.0*1.055) * pow(c, 1.0/2.4) - (255.0*0.055);

	return c;
}


struct cie_XYZ sRGB_to_XYZ(struct RGB sRGB)
{
	double R, G, B;		/* *Linear* RGB values in the tange [0,1] */
	struct cie_XYZ XYZ;

	R = sRGB_to_linear(sRGB.R);
	G = sRGB_to_linear(sRGB.G);
	B = sRGB_to_linear(sRGB.B);

	XYZ.X = 4.12315151515152e-01 * R +
		3.57600000000000e-01 * G +
		1.80500000000000e-01 * B;
	XYZ.Y = 2.12600000000000e-01 * R +
		7.15200000000000e-01 * G +
		7.22000000000000e-02 * B;
	XYZ.Z = 1.93272727272727e-02 * R +
		1.19200000000000e-01 * G +
		9.50633333333333e-01 * B;

	return XYZ;
}

struct RGB const_func XYZ_to_sRGB(struct cie_XYZ XYZ)
{
	double R, G, B;		/* *Linear* RGB values in the tange [0,1] */
	struct RGB sRGB;

	R =   3.24156456493897e+00 * XYZ.X +
	     -1.53766524234284e+00 * XYZ.Y +
	     -4.98702240759841e-01 * XYZ.Z;
	G =  -9.69201189545655e-01 * XYZ.X +
	      1.87588534600745e+00 * XYZ.Y +
	      4.15532375587361e-02 * XYZ.Z;
	B =   5.56241586846000e-02 * XYZ.X +
	     -2.03955248510200e-01 * XYZ.Y +
	      1.05685901500740e+00 * XYZ.Z;

	sRGB.R = linear_to_sRGB(R);
	sRGB.G = linear_to_sRGB(G);
	sRGB.B = linear_to_sRGB(B);

	return sRGB;
}