Calculating RGB plus Amber

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.

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.