Calculating RGB plus Amber

Question

I'm currently working on a wide gamut light source using red, green and blue LED emitters. From an internal xyY (or CIE XYZ) representation, I can reach any color or color temperature via a 3x3 transformation matrix. The matrix is calculated from the chromaticity coordinates and the relative luminance of the actual RGB emitters. This works well both in theory and in practice.

However, the RGB LEDs emitts a discontinuous spectrum with very little energy between red and green. I want to add an orange or amber LED to improve the spectrum and the color rendering index. Starting from CIE xyY, how do I calculate RGB plus Amber? The entire CIE model is based on tristimulus and I can't see how I can use it to calculate a fourth color.

The ideal would probably be a model that could accept any number of colors.

---

ADDED based on discussion:

Above is an illustration of how I imagine the RGB math works.

I measure the chromaticity coordinates (x,y) of each of the Red, Green and Blue emitters, and their relative brightness (Y).

From xyY I calculate CIE XYZ. This is needed because xy plus Y is a projection, XYZ is the actual 3-dimensional color space. I put the XYZ values for the three emitters into a matrix and calculate the inverse of that matrix. This inverse matrix represents the difference between the "actual" (human vision) and the properties of these particular emitters. If I want to display a particular color, say the white shown above to the right, I can take the desired coordinates, send them through the matrix, and get the required strength for each emitter (shown as arrows).

Now RGBA:

Originally I thought that the math for one more emitter (RGBA) was similar to the math for RGB. That I could use a 3*4 matrix to calculate RGBA, the same way I use a 3*3 matrix to calculate RGB. However, it seems like there are problems calculating the inverse of a non-square matrix. Some sources say it can not be done, some say if can be done, but the resulting matrix be lacking in some respect. This is WAY beyond my understanding! :-)

@percusse suggests that a 3*4 matrix can be used. If so, how can I calculate the inverse matrix (I'm on a shoestring budget, software like MATLAB is out of reach).

---

Second EDIT:

Based on the input from @joriki and @percusse I've tried to solve this on paper. I've spent a lot of paper, but I can't seem to find a way to do this that can be implemented as computer code, - or even produces the right answer! I'm probably making mistakes in the practical solving, but that is not actually critical. Computing will have to be done by a C implementation of a solving algorithm (gaussian elimination?) that is known to be good.

This would be typ XYZ values for the four emitters calculated form measured xyY coordinates (actual values will have better precision):

$$ \left[ \begin{array}{ccc} 0.47 & 0.11 & 0.19 & 0.34\\ 0.20 & 0.43 & 0.11 & 0.26\\ 0.00 & 0.06 & 1.12 & 0.01\end{array} \right]

\left[ \begin{array}{ccc} R \\\ G \\\ B \\\ A \end{array} \right]

\left[ \begin{array}{ccc} X \\\ Y \\\ Z \end{array} \right] $$

I've been thinking about optimization and there are a number of parameters that affects the optimal mix, mainly spectrum, efficacy, and heat. For a small system, it is probably enough to worry about the extremes of the amber emitter (avoid max amber when emitting orange-ish light, avoid min amber when emitting any other color). A solution is already suggested by @joriki ["This selects the solution that covers the spectrum most evenly"] but I don't understand the math :-)

So I need to get this system of equations into a form that generates a single answer within 500us of computing time on a small embedded processor :-) Any guidance on how to get a step closer a practical implementation would be greatly appreciated!

---

Third EDIT: I've set up a test that can drive 4 emitters, and a spectrometer to measure the output. The relative intensities of the emitters are tweaked to give a correlated color temperature of roughly 6000 Kelvin (midday daylight).

RGB at ca. 6000K:

RGB + Amber at ca. 6000K:

RGB + White at ca. 6000K:

The first image shows the spectrum from 3 emitters, Red, Green and Blue. There is very little light between 560 and 610 nm. The next image shows the spectrum when Amber is added to RGB. Amber improves the situation significantly. (Yellow might be better, but suitable high brightness yellow LEDs can't be found). The last image shows the spectrum when White is added to RGB. White LEDs are actually phosphor converted blue. The phosphor can be made to retransmitt over a fairly broad spectrum. This seems to give the best result in terms of even spectrum.

I think I have working code for Gaussian Elimination. The question now is how do I add mean square minimization to the equations in such a way that I end up with a single answer? I probably need some hints on how to solve this in practice. Sorry! :-)

---

Fourth and fifth EDIT:

So I have measured the spectra from 380 to 780nm with 1nm resolution. The output is measured at equal input value.

I calculated the area under the curve by trigonometry. I calculated the average size for the 400 trapeziods between 380 and 780nm for R, G, B and A (values are scaled to me more manageable):

$\langle R\rangle = 19.8719507$
$\langle G\rangle = 13.39000051$
$\langle B\rangle = 29.30636046$
$\langle A\rangle = 8.165754589$

And also the average for the product of all six (plus four) combinations of emitter pairs. I then took a stab at assembling this into a covariance matrix:

$$ \left[ \begin{array}{cccc} 43.74282392 & -2.642812728 & -5.823745503 & -0.26554119\\\ -2.642812728 & 8.563382072 & -0.969894212 & -0.946563019\\\ -5.823745503 & -0.969894212 & 62.81754221 & -2.393057209\\\ -0.26554119 & -0.946563019 & -2.393057209 & 8.136438369\end{array} \right] $$

The matrix is assembled like this: $$ \left[ \begin{array}{ccc} \langle RR \rangle - \langle R \rangle\langle R \rangle & \langle RG \rangle - \langle R \rangle\langle G \rangle & \langle RB \rangle - \langle R \rangle\langle B \rangle & \langle RA \rangle - \langle R \rangle\langle A \rangle\\\ \langle RG \rangle - \langle R \rangle\langle G \rangle & ... & ... & ...\\\ ... & ... & ... & ...\\\ ... & ... & ... & ...\end{array} \right] $$

Here are the measured color coordinates of the RGBA emitters, and sample values for XYZ:

$$ \left[ \begin{array}{cccc|c} 0.490449254 & 0.100440581 & 0.221653947 & 0.343906601 & 0.75\\\ 0.204678363 & 0.421052632 & 0.16374269 & 0.210526316 & 1.00\\\ -0.011955512 & 0.07388664 & 1.464803251 & -0.012677086 & 0.75\end{array} \right] $$

I've tried to get the above matrix into echelon form(?) by gaussian elimination, and then get the RGB values on the form $u + Av$ by substitution.

$R: 0.97921341 + A * -0.701207308$
$G: 1.730718699 + A * -0.1767291$
$B: 0.43 + A * 0.012215723$

The next step seems to be to calculate $Q$. This has been answered by @joriki, but am not used to the notation and I'm not at all sure how to translate the greek shorthand to a form where I can calculate the values. If this gets too basic for this forum, let me know and I'll take it offline.

I have trouble understanding this calculation:

$$ \begin{eqnarray} \mu &=& -\frac{\sum_{\alpha,\beta}M_{\alpha\beta}x_\alpha y_\beta}{\sum_{\alpha,\beta}M_{\alpha\beta}y_\alpha y_\beta} \;. \end{eqnarray} $$

Not entirely sure what the $x$ and $y$ values are? A pointer to an example of what this $M_{\alpha\beta}x_\alpha y_\beta$ look like in non-algebraic form would be very helpful.

---

Sixth EDIT:

So let me try to explain how I understand what needs to be done: With a set of measured RGBA emitter color coordinates and an ZYX value (the color we want the emitters to generate) as input we calculate two values for each emitter. The values are

$R = u_{RED} + Av_{RED}$
$G = u_{GREEN} + Av_{GREEN}$
$B = u_{BLUE} + Av_{BLUE}$
$A = A$

The calculation involves gaussian elimination and substitution, and I have written code that performs those calculations.

The value of A should preferably be the one that, together with RGB, produces the most even spectrum. This appears to be the calculation of A:

$$ \begin{eqnarray} A &=& -\frac{\sum_{\alpha,\beta}M_{\alpha\beta}x_\alpha y_\beta}{\sum_{\alpha,\beta}M_{\alpha\beta}y_\alpha y_\beta} \;. \end{eqnarray} $$

One element of this equation is $M_{\alpha\beta}$ which is a 4 * 4 covariance matrix that we have precalculated from the emitter spectra.

This is as far as I am right now. I don't understand from the above notation how the math works. Do I run every possible combination of emitter colors through the matrix add them all up? I have to admit I am completely lost! :-)

Answer 1

As percusse has noted (though using the wrong transformation matrix), the problem is underdetermined. You've now got four variables for three constraints to fulfill.

I presume that when you say you "can reach any colour", you mean that you can specify any tristimulus values within some desired part of the spectrum and calculate the required amplitudes for your R, G and B emitters in order to attain those tristimulus values. Since you're already able to do this without amber, you already know one solution to the underdetermined equations, namely, the one with zero amber.

In order to turn this into a well-defined problem, you need a fourth condition to match the newly introduced fourth variable. Since you're introducing amber in order to cover the spectrum more evenly, one such condition might be to minimize the mean square variation of the intensity over the spectrum.

To do this, you can write the spectrum of your output like this:

$$I(\lambda)=\sum_\alpha c_\alpha I_\alpha(\lambda)\;,$$

where $\alpha$ runs over your four emitters, $c_\alpha$ is the amplitude of emitter $\alpha$, and $I_\alpha(\lambda)$ is the spectrum emitted by emitter $\alpha$. Then the mean square variation is

$$ \def\ave#1{\left\langle#1\right\rangle} \begin{eqnarray} Q &=& \ave{I^2}-\ave{I}^2 \\ &=& \ave{\left(\sum_\alpha c_\alpha I_\alpha\right)^2}-\ave{\sum_\alpha c_\alpha I_\alpha}^2 \\ &=& \sum_{\alpha,\beta}M_{\alpha\beta}c_\alpha c_\beta\;, \end{eqnarray} $$

where

$$\ave f=\frac{\int f(\lambda)\mathrm d\lambda}{\int \mathrm d\lambda}$$

averages over the spectrum and and the covariance matrix $M_{\alpha\beta}$ is given by

$$M_{\alpha\beta}=\ave{I_\alpha I_\beta}-\ave{I_\alpha}\ave{I_\beta}\;.$$

Thus, from the measured spectra of the emitters, you can calculate the covariance matrix, and this gives you a quadratic form that measures how evenly you've covered the spectrum.

You also have a transformation matrix that transforms the emitter amplitudes to tristimulus values, something like

$$d_\gamma=\sum_\alpha T_{\gamma\alpha}c_\alpha\;,$$

where $\gamma$ runs over X, Y, Z and $\alpha$ over the emitters as before. Solving this system of linear equations gives you a one-dimensional solution space $c_\alpha=x_\alpha+\mu y_\alpha$. The positivity constraints on the amplitudes select some interval of admissible $\mu$ values. (Since you say you can already produce all desired colours without amber, this interval isn't empty.) You can substitute this solution into the quality measure $Q$ and minimize:

$$ \begin{eqnarray} Q &=& \sum_{\alpha,\beta}M_{\alpha\beta}c_\alpha c_\beta \\ &=& \sum_{\alpha,\beta}M_{\alpha\beta}(x_\alpha+\mu y_\alpha)(x_\beta+\mu y_\beta)\;, \\ \frac{\mathrm d Q}{\mathrm d\mu} &=& 2\mu\sum_{\alpha,\beta}M_{\alpha\beta}y_\alpha y_\beta + 2\sum_{\alpha,\beta}M_{\alpha\beta}x_\alpha y_\beta \\ &=&0\;, \\ \mu &=& -\frac{\sum_{\alpha,\beta}M_{\alpha\beta}x_\alpha y_\beta}{\sum_{\alpha,\beta}M_{\alpha\beta}y_\alpha y_\beta} \;. \end{eqnarray} $$

This selects the solution that covers the spectrum most evenly, but it may not satisfy the positivity constraints on the amplitudes. If it doesn't, you have to calculate $Q$ at the boundaries of the interval for $\mu$ and use the boundary that yields the lower value.

[Update in response to the edited question:]

It seems my answer assumed too much and thus didn't really answer your question. I'll try to make clearer why this is a solution to your problem.

As percusse has already pointed out in comments, a good way to think about this is in terms of numbers of constraints and numbers of variables. By adding a fourth emitter, you've added a fourth variable: You can now choose what you call the strength of that emitter, in addition to the three RGB strengths you could choose before, and this gives you more choice than you need to get all the colours.

Suppose you have some beer with 4% alcohol and some rum with 40% alcohol. If I ask you to mix me a drink of $X$ ml containing $Y$ ml of alcohol, you can do this by mixing appropriate amounts of the two liquids, as long as $.04X\le Y\le.4 X$ (your "alcohol gamut"), and for given $X$ and $Y$ there is exactly one pair of amounts that you can mix to fulfil that specification. You can determine those amounts by setting up a matrix

$$\pmatrix{1&1\\.04&.4}$$

that transforms from beer content and rum content to total content and alcohol content, and inverting that matrix to transform in the opposite direction.

Now you decide that mixing drinks from beer and rum alone is boring and you get some wine with 12% alcohol. Now your transformation matrix is

$$\pmatrix{1&1&1\\.04&.4&0.12}\;.$$

I think in that case it's clear that there is now no longer a unique solution. For instance, if I ask you for a 100 ml drink with 12 ml of alcohol, you can just give me 100 ml wine, or the mixture of beer and rum that you would have given me before. What's more, you can mix any two solutions and get a new solution; for instance, you can mix 60 ml of wine with 40 ml of the 12% mixture of beer and rum that you would have given me before, and this will also be 100 ml of a 12% drink. This is related to the linearity of the problem: Because the problem is linear, solutions can be linearly superimposed.

Now the fact that the solution is no longer unique is directly related to the fact that a $2\times3$ matrix doesn't have an inverse. If it did, you could use the inverse to find the unique solution. So you can't find the set of solutions by inverting the matrix; but you can still find it, for instance using Gaussian elimination. In the case of the beer, wine and rum, for instance, if I give you a specification of $X$ ml with $Y$ ml of alcohol, you have a linear system of equations:

$$\pmatrix{1&1&1\\.04&.4&0.12}\pmatrix{B\\R\\W}=\pmatrix{X\\Y}\;,$$

where $B$, $R$ and $W$ are the amounts of beer, rum and wine, respectively, to be mixed. Subtracting $.04$ times the upper equation from the lower one yields

$$\pmatrix{1&1&1\\0&.36&0.08}\pmatrix{B\\R\\W}=\pmatrix{X\\Y-.04X}\;.$$

Now you can arbitrarily choose one of the variables, say $W$, and express the others in terms of it:

$$.36R+.08W=Y-.04X\;,$$

$$R=\frac{25}9Y-\frac19X-\frac29W\;,$$

$$B+R+W=X\;,$$

$$B=X-R-W=\frac{10}9X-\frac{25}9Y-\frac79W\;.$$

We can summarize this result as

$$\pmatrix{B\\R\\W}=\pmatrix{\frac{10}9X-\frac{25}9Y\\\frac{25}9Y-\frac19X\\0}+W\pmatrix{-\frac79\\-\frac29\\1}\;.$$

This brings out the linear structure of the solution: You can choose any value of $W$ to get a corresponding solution for $B$, $R$ and $W$ that meets the specification of $X$ and $Y$. However, in general not all of $B$, $R$ and $W$ will be positive, and if they aren't you obviously can't use them as components of a mixture. This is exactly analogous to the triangle in your diagram, which shows the colours you can mix using your emitters; to produce colours outside that triangle, you'd have to be able to run your emitters with negative strengths.

So to summarize, adding a fourth emitter has given you an extra degree of freedom; there's no longer a unique solution but a whole family of solutions that depends linearly on some parameter. Thus, you need a criterion to pick out one of these answers.

The thinking behind my answer was: This criterion should have something to do with why you decided to introduce the amber emitter in the first place – otherwise you might as well stick with the solution you had before, which is of course still a solution since you can always choose strength $0$ for the amber emitter. Since you said you added the amber emitter because the existing emitters were producing a discontinuous spectrum, I figured that a useful criterion might be to minimize the variance in the spectrum.

One more point: In your edit you only talk about the XYZ measurements, which give you the transformation matrix. But I presume you must have also made full spectrum measurements, since otherwise you wouldn't know that the spectrum isn't smooth. You're going to need those measurements to apply the solution I proposed, or in fact any solution that uses a criterion related to the smoothness of the spectrum to fix the extra parameter.

[Second edit]

Your description of how you'd like to pick one of the possible solutions isn't really concrete enough to offer a concrete solution, but as a first step here's how to find the solution set in your case by Gaussian elimination. I'll stick with two digits precision (though it leads to significant rounding errors in this case).

$$ \pmatrix{ 0.47 & 0.11 & 0.19 & 0.34\\ 0.20 & 0.43 & 0.11 & 0.26\\ 0.00 & 0.06 & 1.12 & 0.01 } \pmatrix{ R \\ G \\ B \\ A } = \pmatrix{ X \\ Y \\ Z } $$

Subtract $0.20/0.47$ times the first row from the second:

$$ \pmatrix{ 0.47 & 0.11 & 0.19 & 0.34\\ 0.00 & 0.38 & 0.03 & 0.12\\ 0.00 & 0.06 & 1.12 & 0.01 } \pmatrix{ R \\ G \\ B \\ A } = \pmatrix{ X \\ Y-0.43X \\ Z } $$

Subtract $0.06/0.38$ times the second row from the third:

$$ \pmatrix{ 0.47 & 0.11 & 0.19 & 0.34\\ 0.00 & 0.38 & 0.03 & 0.12\\ 0.00 & 0.00 & 1.12 & -0.01 } \pmatrix{ R \\ G \\ B \\ A } = \pmatrix{ X \\ Y-0.43X \\ Z-0.07X-0.16Y } $$

Now we can freely choose $A$ and express the other variables in terms of $A$. Solve the third row for $B$:

$$1.12B-0.01A=Z-0.07X-0.16Y\;,$$ $$B=\frac{Z-0.07X-0.16Y+0.01A}{1.12}\;.$$

Analogously solve the second row for $G$ and the first row for $R$. Then collect the terms independent of $A$ and the terms proportional to $A$ in separate vectors:

$$\pmatrix{R\\G\\B\\A}=\pmatrix{\ldots\\\ldots\\(Z-0.07X-0.16Y)/1.12\\0}+A\pmatrix{\ldots\\\ldots\\0.01/1.12\\1}\;,$$

where the ellipses stand for results you obtain from solving the first and second rows. For given $X$, $Y$ and $Z$, this is of the form $u+Av$ with constant vectors $u$ and $v$. You can find the interval of admissible $A$ values that make all emitter strengths non-negative, and then you can apply whatever criterion you decide on to select $A$ within that interval.

[Third edit]

Hi Morten -- your edit made me realize I forgot to respond to your email with the code -- sorry about that! Here's what you need to do now:

To calculate the matrix $M_{\alpha\beta}$, you need to calculate integrals over the spectrum. $\ave{I_\alpha}$ is the average over the spectrum for emitter $\alpha$, and $\ave{I_\alpha I_\beta}$ is the average over the product of the spectra for emitters $\alpha$ and $\beta$. When you've calculated these, you can assemble $M_{\alpha\beta}$ from them, and then together with the vectors $x_\alpha$ and $y_\alpha$ you get out of the Gaussian elimination you have all the ingredients for the formula for the optimal parameter $\mu$. Once you have that, you need to test whether it yields non-negative emitter strengths. If it does, you're done; if it doesn't, you need to find the boundaries of the interval of values that do, and pick the one that leads to the lower $Q$ value.

I'm aware that there may be details in some of those steps that you'll need further help with, but I'll leave it at that sketch for now and you can say specifically where you'd like me to fill in more details.

Answer 2

If I understand correctly, your problem can be modeled as the following:

There is a linear bijective transformation between CIE XYZ space and RGB space via (stealing from this wiki page, heh) $$ \begin{bmatrix}X\\Y\\Z\end{bmatrix}=\frac{1}{0.17697}\begin{bmatrix}0.49&0.31&0.20\\0.17697&0.81240&0.01063\\0.00&0.01&0.99\end{bmatrix}\begin{bmatrix}R\\G\\B\end{bmatrix} $$ This is the common $Ar=x$ case with an invertible $A$. The inverse is simply $A^{-1}x=r$

What you want to achieve is,to solve the following equation.

Given $XYZ$ data and the mapping $M$, compute me the 4 entry vector including RGBAm values. Writing it down, $$ \underbrace{\pmatrix{\cdot&\cdot&\cdot&a\\\cdot&\cdot&\cdot&b\\\cdot&\cdot&\cdot&c}}_{M}\begin{bmatrix}R\\G\\B\\Am\end{bmatrix} = \begin{bmatrix}X\\Y\\Z\end{bmatrix} $$

Here we have the $Am$ blues, in other words, you must provide additional information. The missing piece is the numerical values of $a,b,c$. In physical terms you have to quantify the effect of amber at each $X,Y,Z$.

Now, my tiny knowledge of light tells me that Amber must be a combination of RGB anyway. So the amber effect can be quantized as $Am = aR+bG+cB$. And this can info can be obtained by the Amber led specifications you are using.

The remaining step is the well-known intensely studied and also well-established underdetermined least squares problem. Note that now, we have infinitely many (well,theoretically that is) $R,G,B,Am$ quartets that would lead to the same $XYZ$. In other words the map is not one-to-one (bijective) any more but rather onto (surjective). The solution is the MATLAB code linsolve(M,x) for you to tryout.

Assume this as a proof of concept type of answer. Please correct my misunderstanding and feel free to edit until we arrive to a concrete problem. Also my sincere congratulations for your project. It is a fantastic initiative.