summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
authorH. Peter Anvin <hpa@zytor.com>2017-09-12 18:11:41 (GMT)
committerH. Peter Anvin <hpa@zytor.com>2017-09-12 18:26:50 (GMT)
commit7b31610abcf59c85b31a0b7c3db3a33ed82f71fd (patch)
treead7fac9c532e798a3dd178994d784d91595d9765
parentb92ea995371cec99e0366f6a341ca1e35e1f7211 (diff)
downloadledcolor-master.zip
ledcolor-master.tar.gz
ledcolor-master.tar.bz2
ledcolor-master.tar.xz
Initialize RGB <-> XYZ matrices in C code, to make other LEDs easierHEADmaster
Initialize the RGB <-> XYZ matrices explicitly in C code, eliminating the need to use octave, and making it easier to specify different primaries for other LEDs.
-rw-r--r--Makefile2
-rw-r--r--color.c34
-rw-r--r--color.h16
-rw-r--r--matrix.c57
-rw-r--r--mktables.c10
-rw-r--r--rgbxyz.c107
-rw-r--r--rgbxyz.mat20
-rw-r--r--spectral.csv42
8 files changed, 213 insertions, 75 deletions
diff --git a/Makefile b/Makefile
index 069393e..ddf28b0 100644
--- a/Makefile
+++ b/Makefile
@@ -11,7 +11,7 @@ csv: mktables
%.o: %.c
$(CC) $(CFLAGS) -c -o $@ $<
-mktables: mktables.o color.o rgbxyz.o
+mktables: mktables.o color.o rgbxyz.o matrix.o
$(CC) $(LDLFLAGS) -o $@ $^ $(LIBS)
clean:
diff --git a/color.c b/color.c
index f50d51f..ce34b5b 100644
--- a/color.c
+++ b/color.c
@@ -115,18 +115,42 @@ struct cie_XYZ const_func Luv_to_XYZ(struct cie_Luv Luv, struct cie_XYZ white)
struct cie_xy pure_func spectral(double nm)
{
- /* If the value is black, return the point of equal power */
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];
+ 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)
*/
@@ -191,13 +215,13 @@ static double const_func linear_to_sRGB(double c)
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);
@@ -219,7 +243,7 @@ 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;
diff --git a/color.h b/color.h
index e4811c4..7c7ab8b 100644
--- a/color.h
+++ b/color.h
@@ -54,9 +54,19 @@ struct cie_Luv const_func XYZ_to_Luv(struct cie_XYZ XYZ, struct cie_XYZ white);
struct cie_XYZ const_func Luv_to_XYZ(struct cie_Luv Luv, struct cie_XYZ white);
struct cie_xy pure_func spectral(double nm);
struct cie_XYZ pure_func planckian(double T);
-struct RGB const_func XYZ_to_RGB(struct cie_XYZ XYZ);
-struct cie_XYZ const_func RGB_to_XYZ(struct RGB RGB);
+struct RGB pure_func XYZ_to_RGB(struct cie_XYZ XYZ);
+struct cie_XYZ pure_func RGB_to_XYZ(struct RGB RGB);
+void RGBXYZ_init(void);
struct RGB const_func XYZ_to_sRGB(struct cie_XYZ XYZ);
-struct cie_XYZ sRGB_to_XYZ(struct RGB sRGB);
+struct cie_XYZ const_func sRGB_to_XYZ(struct RGB sRGB);
+/* y = A x */
+static inline void color_mat_multiply(double y[3], const double A[3][3],
+ const double x[3])
+{
+ y[0] = A[0][0] * x[0] + A[0][1] * x[1] + A[0][2] * x[2];
+ y[1] = A[1][0] * x[0] + A[1][1] * x[1] + A[1][2] * x[2];
+ y[2] = A[2][0] * x[0] + A[2][1] * x[1] + A[2][2] * x[2];
+}
+void color_mat_inverse(double X[3][3], const double A[3][3]); /* X = A^-1 */
#endif /* COLOR_H */
diff --git a/matrix.c b/matrix.c
new file mode 100644
index 0000000..3d86681
--- /dev/null
+++ b/matrix.c
@@ -0,0 +1,57 @@
+/*
+ * Matrix inversion for a 3x3 matrix
+ */
+
+#include "color.h"
+
+/*
+ * X = A^-1
+ *
+ * This uses the "method of adjugate" to calculate the inverse
+ * matrix. If the matrix is singular (meaning the columns are
+ * linearly equivalent "in a straight line") then the result will
+ * contain infinities.
+ */
+void color_mat_inverse(double X[3][3], const double A[3][3])
+{
+ double a[3][3]; /* Adjugate matrix */
+ double det; /* Determinant */
+ double invdet; /* Inverse determinant */
+ int i, j;
+
+ /*
+ * The "matrix of cofactors" is the determinate of each
+ * individual submatrix ignoring the row and column for
+ * which the value is computed, multiplied with (-1)^(i+j)
+ *
+ * The "adjugate" is the transpose of this matrix.
+ */
+
+ a[0][0] = A[1][1]*A[2][2] - A[1][2]*A[2][1]; /* \ - / */
+ a[0][1] = A[0][2]*A[2][1] - A[0][1]*A[2][2]; /* / - \ */
+ a[0][2] = A[0][1]*A[1][2] - A[0][2]*A[1][1]; /* \ - / */
+ a[1][0] = A[1][2]*A[2][0] - A[1][0]*A[2][2]; /* / - \ */
+ a[1][1] = A[0][0]*A[2][2] - A[0][2]*A[2][0]; /* \ - / */
+ a[1][2] = A[0][2]*A[1][0] - A[0][0]*A[1][2]; /* / - \ */
+ a[2][0] = A[1][0]*A[2][1] - A[1][1]*A[2][0]; /* \ - / */
+ a[2][1] = A[0][1]*A[2][0] - A[0][0]*A[2][1]; /* / - \ */
+ a[2][2] = A[0][0]*A[1][1] - A[0][1]*A[1][0]; /* \ - / */
+
+ /*
+ * The determinant is now simply the inner product of the
+ * first column of the adjugate with top row of the original
+ * matrix (or any other row/column)
+ */
+ det = A[0][0]*a[0][0] + A[1][0]*a[0][1] + A[2][0]*a[0][2];
+ invdet = 1.0/det;
+
+ /*
+ * The inverse is the adjugate divided with the determinant.
+ * If the matrix is singular (non-invertible), then
+ * det = 0 and invdet = infinity, so we get a matrix of
+ * infinities here.
+ */
+ for (i = 0; i < 3; i++)
+ for (j = 0; j < 3; j++)
+ X[i][j] = a[i][j] * invdet;
+}
diff --git a/mktables.c b/mktables.c
index f5ac4cd..015fd99 100644
--- a/mktables.c
+++ b/mktables.c
@@ -183,7 +183,7 @@ static void spectral_table(FILE *f)
RGB = XYZ_to_RGB(XYZ);
printf("Peak: RGB = %.1f %.1f %.1f Y = %.6f xy = %.3f %.3f\n",
RGB.R, RGB.G, RGB.B, XYZ.Y, xy.x, xy.y);
-
+
RGB.R = 255.0;
RGB.G = 0.0;
RGB.B = 0.0;
@@ -192,7 +192,7 @@ static void spectral_table(FILE *f)
RGB = XYZ_to_RGB(XYZ);
printf("Redk: RGB = %.1f %.1f %.1f Y = %.6f xy = %.3f %.3f\n",
RGB.R, RGB.G, RGB.B, XYZ.Y, redxy.x, redxy.y);
-
+
for (nm = 466.0; nm <= 800.05; nm += 0.1) {
xy = spectral(nm);
XYZ = xyY_to_XYZ(xy, blue.Y);
@@ -218,7 +218,11 @@ static void spectral_table(FILE *f)
*/
int main(void)
{
- FILE *f = fopen("thermal.csv", "w");
+ FILE *f;
+
+ RGBXYZ_init();
+
+ f = fopen("thermal.csv", "w");
thermal_table(f);
fclose(f);
diff --git a/rgbxyz.c b/rgbxyz.c
index af8bdc9..c8de1a0 100644
--- a/rgbxyz.c
+++ b/rgbxyz.c
@@ -1,10 +1,22 @@
/*
- * Matrixes for conversion between CIE XYZ values to WS2182B RGB,
+ * Matrixes for conversion between CIE XYZ values to device linear RGB,
* with values normalized so that Y = 1.0 and R = G = B = 255.0
* represent the maximum possible values.
+ *
+ * If the device primaries are nonlinear, or if there are more then
+ * three primaries, additional code could be added to these functions,
+ * possibly replacing the type of "struct RGB". Either way, such code
+ * would involve more than a simple matrix multiplication;
*/
#include "color.h"
+#include <math.h>
+#include <stdio.h>
+
+/*
+ * Conversion matrices, calculated by RGBXYZ_init()
+ */
+static double RGBXYZ[3][3], XYZRGB[3][3];
/*
* This is the inverse of the RGB-to-XYZ matrix
@@ -13,18 +25,7 @@ struct RGB const_func XYZ_to_RGB(struct cie_XYZ XYZ)
{
struct RGB RGB;
- RGB.R = ( 3.96394821253806e+02) * XYZ.X +
- (-7.43756319424205e+01) * XYZ.Y +
- (-5.39829972623372e+01) * XYZ.Z;
-
- RGB.G = (-2.33573501299311e+02) * XYZ.X +
- ( 5.20567953949116e+02) * XYZ.Y +
- (-2.29114388822662e+01) * XYZ.Z;
-
- RGB.B = ( 3.29375415122602e+00) * XYZ.X +
- (-7.34082783268273e+00) * XYZ.Y +
- ( 1.98273719579001e+02) * XYZ.Z;
-
+ color_mat_multiply((double *)&RGB, XYZRGB, (const double *)&XYZ);
return RGB;
}
@@ -45,16 +46,78 @@ struct cie_XYZ const_func RGB_to_XYZ(struct RGB RGB)
{
struct cie_XYZ XYZ;
- XYZ.X = ( 2.75464459312007e-03) * RGB.R +
- ( 4.04802838039216e-04) * RGB.G +
- ( 7.96770178823529e-04) * RGB.B;
+ color_mat_multiply((double *)&XYZ, RGBXYZ, (const double *)&RGB);
+ return XYZ;
+}
- XYZ.Y = ( 1.23598077372549e-03) * RGB.R +
- ( 2.10574502156863e-03) * RGB.G +
- ( 5.79842832156863e-04) * RGB.B;
+/*
+ * Compute the RGB <-> XYZ transformation matrices. For pure spectral
+ * primaries (as one can expect from LEDs or lasers) this code should
+ * be possible to use as-is with modified constants; for nonspectral
+ * primaries the XYZ or xyY values will have to be obtained directly
+ * some other way.
+ */
+void RGBXYZ_init(void)
+{
+ /*
+ * This code assumes the available data specifies a luminous
+ * (photometric) intensity, normally given in candela (cd).
+ * If what we have is a bolometric intensity (W/sr) then
+ * the values need to be scaled using the photometric_factor()
+ * for the given wavelength.
+ */
+ struct primary {
+ double nm; /* Wavelength in nm */
+ double cd; /* Intensity in cd */
+ };
- XYZ.Z = ( 7.12378356862745e-05) * RGB.G +
- ( 5.05176461843137e-03) * RGB.B;
+ /* From the WS2812B data sheet (average of ranges given) */
+ const struct primary primary[3] = {
+ { 622.5, 0.405 }, /* Red */
+ { 523.5, 0.690 }, /* Green */
+ { 466.0, 0.190 } /* Blue */
+ };
+ double peakcd; /* Peak (total) intensity in cd */
+ double scale;
+ int i;
- return XYZ;
+ peakcd = 0.0;
+ for (i = 0; i < 3; i++)
+ peakcd += primary[i].cd;
+
+ /*
+ * Compute the RGBXYZ matrix. Note that we already have a
+ * photometric intensity (in cd). If we have a *bolometric*
+ * intensity (in W/sr) then we need to convert it to a photometric
+ * intensity using the photometric_scale() conversion factor.
+ *
+ * This matrix represents the contribution to the total XYZ
+ * value for *each scale point* of each primary.
+ */
+ scale = 1.0/(255.0*peakcd);
+
+ printf("rgbxyz =\n\n");
+ for (i = 0; i < 3; i++) {
+ struct cie_xy xy = spectral(primary[i].nm);
+ struct cie_XYZ XYZ = xyY_to_XYZ(xy, primary[i].cd*scale);
+
+ RGBXYZ[0][i] = XYZ.X;
+ RGBXYZ[1][i] = XYZ.Y;
+ RGBXYZ[2][i] = XYZ.Z;
+ }
+ for (i = 0; i < 3; i++) {
+ printf(" %22.14e %22.14e %22.14e\n",
+ RGBXYZ[i][0], RGBXYZ[i][1], RGBXYZ[i][2]);
+ }
+
+ /*
+ * The inverse matrix provides for XYZ -> RGB conversion
+ */
+ color_mat_inverse(XYZRGB, RGBXYZ);
+ printf("\nxyzrgb =\n\n");
+ for (i = 0; i < 3; i++) {
+ printf(" %22.14e %22.14e %22.14e\n",
+ XYZRGB[i][0], XYZRGB[i][1], XYZRGB[i][2]);
+ }
+ putchar('\n');
}
diff --git a/rgbxyz.mat b/rgbxyz.mat
deleted file mode 100644
index e5e6235..0000000
--- a/rgbxyz.mat
+++ /dev/null
@@ -1,20 +0,0 @@
-# Octave/Matlab script for computing the RGB <-> XYZ matrices
-
-# Values from CIE table @ 622.5 nm, 523.5 nm, 466.0 nm
-rgbxyz = transpose([[8.134280E-01,3.649768E-01,0.000000E+00],
- [1.580186E-01,8.219974E-01,2.780836E-02],
- [2.420135E-01,1.761233E-01,1.534439E+00]])
-
-# Scale to 405 mcd, 690 mcd, 190 mcd normalized to Y = 1 at
-# total = 1285 mcd
-rgbxyz(1:3,1) *= 405/(1285*255) / rgbxyz(2,1)
-rgbxyz(1:3,2) *= 690/(1285*255) / rgbxyz(2,2)
-rgbxyz(1:3,3) *= 190/(1285*255) / rgbxyz(2,3)
-
-xyzrgb = inverse(rgbxyz)
-
-format long e
-rgbxyz
-xyzrgb
-
-
diff --git a/spectral.csv b/spectral.csv
index fcac437..68876f4 100644
--- a/spectral.csv
+++ b/spectral.csv
@@ -27,7 +27,7 @@
468.6,0.000000,14.555498,212.589794
468.7,0.000000,15.083526,211.040829
468.8,0.000000,15.609033,209.498245
-468.9,0.000000,16.132015,207.962266
+468.9,0.000000,16.132015,207.962265
469.0,0.000000,16.652445,206.432749
469.1,0.000000,17.170313,204.909789
469.2,0.000000,17.685613,203.393419
@@ -104,7 +104,7 @@
476.3,0.000000,47.451463,113.583854
476.4,0.000000,47.775752,112.575746
476.5,0.000000,48.097375,111.575053
-476.6,0.000000,48.416362,110.581710
+476.6,0.000000,48.416362,110.581709
476.7,0.000000,48.732680,109.595839
476.8,0.000000,49.046337,108.617389
476.9,0.000000,49.357340,107.646374
@@ -146,7 +146,7 @@
480.5,0.000000,58.820944,77.585360
480.6,0.000000,59.037907,76.881327
480.7,0.000000,59.252618,76.183737
-480.8,0.000000,59.465091,75.492548
+480.8,0.000000,59.465091,75.492547
480.9,0.000000,59.675367,74.807612
481.0,0.000000,59.883463,74.128899
481.1,0.000000,60.089410,73.456291
@@ -159,7 +159,7 @@
481.8,0.000000,61.472763,68.913375
481.9,0.000000,61.662340,68.287111
482.0,0.000000,61.849968,67.666338
-482.1,0.000000,62.035670,67.050977
+482.1,0.000000,62.035669,67.050977
482.2,0.000000,62.219470,66.440938
482.3,0.000000,62.401383,65.836198
482.4,0.000000,62.581428,65.236679
@@ -238,7 +238,7 @@
489.7,0.000000,71.593883,32.689229
489.8,0.000000,71.669663,32.372791
489.9,0.000000,71.744377,32.059336
-490.0,0.000000,71.818033,31.748842
+490.0,0.000000,71.818032,31.748842
490.1,0.000000,71.890642,31.441276
490.2,0.000000,71.962212,31.136620
490.3,0.000000,72.032749,30.834850
@@ -299,7 +299,7 @@
495.8,0.000000,74.520654,17.976702
495.9,0.000000,74.543830,17.799479
496.0,0.000000,74.566382,17.623844
-496.1,0.000000,74.588323,17.449774
+496.1,0.000000,74.588322,17.449774
496.2,0.000000,74.609665,17.277240
496.3,0.000000,74.630421,17.106227
496.4,0.000000,74.650604,16.936700
@@ -345,11 +345,11 @@
500.4,0.000000,75.115135,11.162356
500.5,0.000000,75.120868,11.038990
500.6,0.000000,75.126390,10.916515
-500.7,0.000000,75.131697,10.794938
+500.7,0.000000,75.131696,10.794938
500.8,0.000000,75.136787,10.674251
500.9,0.000000,75.141659,10.554454
501.0,0.000000,75.146311,10.435543
-501.1,0.000000,75.150740,10.317513
+501.1,0.000000,75.150739,10.317513
501.2,0.000000,75.154943,10.200369
501.3,0.000000,75.158920,10.084103
501.4,0.000000,75.162668,9.968714
@@ -1004,7 +1004,7 @@
566.3,38.345558,47.892263,0.000000
566.4,38.474681,47.816215,0.000000
566.5,38.604189,47.739938,0.000000
-566.6,38.734090,47.663431,0.000000
+566.6,38.734089,47.663431,0.000000
566.7,38.864362,47.586704,0.000000
566.8,38.994996,47.509763,0.000000
566.9,39.125985,47.432613,0.000000
@@ -1012,7 +1012,7 @@
567.1,39.388993,47.277706,0.000000
567.2,39.520989,47.199962,0.000000
567.3,39.653301,47.122030,0.000000
-567.4,39.785926,47.043915,0.000000
+567.4,39.785926,47.043914,0.000000
567.5,39.918839,46.965629,0.000000
567.6,40.052048,46.887169,0.000000
567.7,40.185536,46.808543,0.000000
@@ -1027,7 +1027,7 @@
568.6,41.397947,46.094410,0.000000
568.7,41.533721,46.014435,0.000000
568.8,41.669677,45.934352,0.000000
-568.9,41.805794,45.854175,0.000000
+568.9,41.805793,45.854175,0.000000
569.0,41.942079,45.773897,0.000000
569.1,42.078512,45.693533,0.000000
569.2,42.215089,45.613083,0.000000
@@ -1135,7 +1135,7 @@
579.4,57.104032,36.841612,0.000000
579.5,57.266302,36.746005,0.000000
579.6,57.429113,36.650080,0.000000
-579.7,57.592275,36.553949,0.000000
+579.7,57.592274,36.553948,0.000000
579.8,57.755974,36.457500,0.000000
579.9,57.920093,36.360803,0.000000
580.0,58.084672,36.263836,0.000000
@@ -1300,7 +1300,7 @@
595.9,83.977076,21.007724,0.000000
596.0,84.137335,20.913295,0.000000
596.1,84.297475,20.818935,0.000000
-596.2,84.457471,20.724662,0.000000
+596.2,84.457471,20.724661,0.000000
596.3,84.617499,20.630368,0.000000
596.4,84.777371,20.536167,0.000000
596.5,84.937277,20.441945,0.000000
@@ -1359,7 +1359,7 @@
601.8,93.390681,15.460922,0.000000
601.9,93.549498,15.367341,0.000000
602.0,93.708087,15.273895,0.000000
-602.1,93.866538,15.180531,0.000000
+602.1,93.866537,15.180531,0.000000
602.2,94.024867,15.087237,0.000000
602.3,94.182967,14.994079,0.000000
602.4,94.340885,14.901028,0.000000
@@ -1421,7 +1421,7 @@
608.0,102.652676,10.003415,0.000000
608.1,102.790347,9.922294,0.000000
608.2,102.927561,9.841442,0.000000
-608.3,103.064550,9.760723,0.000000
+608.3,103.064549,9.760723,0.000000
608.4,103.201076,9.680277,0.000000
608.5,103.337244,9.600041,0.000000
608.6,103.473138,9.519968,0.000000
@@ -1465,7 +1465,7 @@
612.4,108.394075,6.620361,0.000000
612.5,108.517478,6.547647,0.000000
612.6,108.640561,6.475122,0.000000
-612.7,108.763320,6.402788,0.000000
+612.7,108.763319,6.402788,0.000000
612.8,108.885776,6.330631,0.000000
612.9,109.007903,6.258670,0.000000
613.0,109.129805,6.186840,0.000000
@@ -1482,7 +1482,7 @@
614.1,110.450511,5.408628,0.000000
614.2,110.568776,5.338942,0.000000
614.3,110.686753,5.269425,0.000000
-614.4,110.804442,5.200078,0.000000
+614.4,110.804441,5.200078,0.000000
614.5,110.921810,5.130920,0.000000
614.6,111.038896,5.061928,0.000000
614.7,111.155695,4.993105,0.000000
@@ -1505,7 +1505,7 @@
616.4,113.098769,3.848274,0.000000
616.5,113.210835,3.782240,0.000000
616.6,113.322684,3.716334,0.000000
-616.7,113.434325,3.650550,0.000000
+616.7,113.434324,3.650550,0.000000
616.8,113.545752,3.584892,0.000000
616.9,113.656961,3.519363,0.000000
617.0,113.767991,3.453939,0.000000
@@ -1533,18 +1533,18 @@
619.2,116.168297,2.039572,0.000000
619.3,116.275848,1.976198,0.000000
619.4,116.383277,1.912896,0.000000
-619.5,116.490625,1.849642,0.000000
+619.5,116.490624,1.849642,0.000000
619.6,116.597856,1.786456,0.000000
619.7,116.704978,1.723335,0.000000
619.8,116.812021,1.660261,0.000000
619.9,116.918992,1.597228,0.000000
-620.0,117.025842,1.534268,0.000000
+620.0,117.025841,1.534268,0.000000
620.1,117.132627,1.471345,0.000000
620.2,117.239268,1.408508,0.000000
620.3,117.345832,1.345715,0.000000
620.4,117.452255,1.283006,0.000000
620.5,117.558538,1.220379,0.000000
-620.6,117.664653,1.157852,0.000000
+620.6,117.664652,1.157852,0.000000
620.7,117.770565,1.095444,0.000000
620.8,117.876310,1.033134,0.000000
620.9,117.981805,0.970972,0.000000