Color: Cohen’s intriguing F matrix again

Let me run thru the derivation of Figure 20 again, still with the CIE 1964 tables, but at 20 nm intervals. Then I show two alternative calculations. And one of them will show us why Cohen switched the sign on f2.

There is interesting linear algebra and vector algebra in this post, even if you’re not interested in where these problems came from.

Review:

What we are aiming for is an orthonormal basis whose three vectors are the columns of an F matrix. In particular, the first column — the first basis vector — is the fundamental of an equal-energy illuminant. The challenge was to find two other columns of F which matched Cohen’s Figure 20.

There were scaling issues all over the place. In particular, Cohen’s final three vectors were orthogonal but not orthonormal, because they were not of unit length. Having matched Figure 20, I will no longer concern myself with his scaling; unit vectors, I want unit vectors.

Let’s grab things as we need them.

Since I’ve shown you pictures derived from the 1964 xyz bar tables, I’ll stay with 1964.

What range of wavelengths? The CIE 1964 xyz bar tables run from 380 to 780. What the heck, just use 400 to 700 as before.

What interval? 20 nm, so I can display things.

Then my A matrix is the following subset of the CIE 1964 xyz bar tables:

There are 16 distinct wavelengths. I take the equal-energy illuminant to have the value 100 for 16 values:

ee= {100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100}

To get its fundamental, I need a projection operator onto the fundamental color space (the pre-image of the range of A’, the transpose of A). To get that, I need an orthonormal basis for the fundamental color space. To get that, I need the Singular Value Decomposition (SVD) of A’:

A’ = u w v’

We know that the first 3 columns (because A only had 3 columns) of v are an orthonormal basis for the fundamental color space; the remaining 13 (= 16 – 3) columns of v are an orthonormal basis for the null space of A’. Here are the first 3 columns of v, and I name them v1:

We know that we construct a projection operator R onto the column space of v1 by

R = v1 v1′.

We can quickly check that R is, at least, a projection operator, since it is idempotent: i.e. R^2 = R.

And we can check that it is a projection operator onto the column space of v1, since it preserves each of the columns of v1: i.e. R v1 = v1.

If we wanted, we could confirm that it projects all of the other columns of v to the zero vector: i.e. R v2 = 0, where v2 is the last 13 columns of v.

Enough already, we came here to get the fundamental of the equal energy illuminant…

f = R ee = {5.53946, 62.01, 125.887, 116.934, 68.0773, 55.4377, 80.9502, 109.171, 126.547, 128.189, 113.072, 76.3005, 36.473, 12.5844, 3.34797, 0.782783}

Make it a unit vector f1…

f1 = {0.0165427, 0.185182, 0.375939, 0.349205, 0.203302, 0.165555, 0.241744, 0.326021, 0.377912, 0.382816, 0.337672, 0.227859, 0.108921, 0.0375811, 0.00999815, 0.00233765}

For the other two vectors, we found that we could match Cohen by taking the three vectors {f1, E2, E3} and applying the Gram-Schmidt orthogonalization procedure… which is nothing more than subtracting off the part of E2 which is parallel to f1 (and normalizing, to get a unit vector f2), and then subtracting off the part of E3 which is in the span of f1 and f2.

First, we need the E basis; they are the dual basis to A’; or E’ is the dual basis to A. As always, I find the solution by writing the definition…

E’A = I

then writing the false (because A’ is not invertible)…

,

which leads me to write the true

.

As usual, that formula is for the transpose E’ but I write E = E”, the transpose of the transpose:

Just for my sanity, let me make unit vectors e2 and e3 out of E2 and E3 (which are, in case you have any doubt, the 2nd and 3rd columns of E):

e2 = {-0.00645091, -0.0664311, -0.0997151, -0.00489041, 0.190674, 0.378755, 0.53808, 0.5274, 0.34079, 0.0360794, -0.204679, -0.238978, -0.137198, -0.0511362, -0.0139157, -0.003292}

e3 = {0.02679, 0.305381, 0.626647, 0.572896, 0.28794, 0.130352, 0.0827338, 0.0307357, -0.0437682, -0.12778, -0.176551, -0.146361, -0.0761614, -0.0273034, -0.00734788, -0.00172836}

The part of e2 which is parallel to f1 is…

f1 (e2.f1)
= {0.00585458, 0.0655375, 0.133048, 0.123587, 0.0719501, 0.0585914, 0.0855552, 0.115382, 0.133746, 0.135482, 0.119505, 0.080641, 0.0385479, 0.0133002, 0.00353843, 0.000827313}

so the part of e2 which is orthogonal to f1 is found by subtracting the parallel part from e2…

t1 = e2-f1 (e2.f1)
= {-0.0123055, -0.131969, -0.232763, -0.128477, 0.118724, 0.320164, 0.452525, 0.412018, 0.207044, -0.0994022, -0.324183, -0.319619, -0.175746, -0.0644364, -0.0174541, -0.00411931}

Do a quick check that t1 is othogonal to f1…

t1.f1 = -5.15775*10^-17

But t1 is not a unit vector: it’s squared length is t1.t1 = 0.874749 .

So we normalize it, and call the result f2:

f2 = {-0.013157, -0.141101, -0.24887, -0.137367, 0.126939, 0.342318, 0.483839, 0.440529, 0.221371, -0.106281, -0.346616, -0.341736, -0.187908, -0.0688953, -0.0186619, -0.00440436}

To get f3, we find the part of e3 which lies in the plane spanned by f1 and f2…

f1 (f1.e3)+f2 (f2.e3)
= {0.00735195, 0.0822258, 0.166479, 0.153523, 0.0869147, 0.0679319, 0.0993851, 0.136557, 0.161737, 0.167772, 0.150999, 0.103177, 0.0496131, 0.0171664, 0.00457096, 0.00106922}

then subtract that from e3, and the result will be orthogonal to the plane spanned by f1 and f2…

t2 = e3-(f1 (f1.e3)+f2 (f2.e3))
= {0.0194381, 0.223156, 0.460167, 0.419372, 0.201026, 0.0624205, -0.0166512, -0.105822, -0.205505, -0.295552, -0.32755, -0.249538, -0.125774, -0.0444698, -0.0119188, -0.00279758}

Like t1, t2 is not a unit vector, but we make it so, and call the result f3:

f3 = {0.0215889, 0.247847, 0.511084, 0.465775, 0.223269, 0.0693272, -0.0184937, -0.117531, -0.228244, -0.328255, -0.363793, -0.277149, -0.139691, -0.0493903, -0.0132376, -0.00310713}

If we arrange f1, f2, f3 as the columns of F…

we can compute F’F. Since F is supposed to be orthonormal, we expect to get a 3×3 identity matrix. Sure enough…

The resulting F, of course, are far from smooth…

but they still have the same kind of shapes we saw in Figure 20.

Oh, except that to match his figure, I had to take the negative of f2.

But I know of no good reason to do that, so I’m not going to. Yet. (Ok, Ok — the obvious guess is that f1, f2, f3 is a left-handed basis. If we change the sign of f2, we will get a right-handed basis. I’ll show it to you later in this post. That is, there is a good reason to do that, but nothing we’ve seen yet would lead us to it. )

What have we done? On a coarser scale, we have reproduced the calculations that lead to a specific orthonormal F matrix. It is one of many — infinitely many — bases which have f1 as the first basis vector. That really is the only special vector; the other two span the plane perpendicular to f1.

Cohen’s Cholesky Decomposition

Now let’s try it the way I think Cohen did it.

There is one potentially confusing problem: notation, as so often happens.

In the following discussion, “A” is no longer a matrix of color-matching functions. We will take “A” to have columns f, E2, E3; that was my starting point for the Gram-Schmidt orthogonalization. He describes such a starting point, but I believe he never shows the computations for such a case.

Here then is our new A matrix (which I call A1 in practice):

Cohen computed two F matrices (F1 and F2) using two Cholesky decompositions, but we only need one of them. We recall that the first column of F2 = E Ga is proportional to the first column of A — and that’s what we want — so we will need E1, i.e. the dual basis to A1.

That is crucial. We need the dual basis to A1, not to A. We’ve already used the dual basis to A, to get two columns of A1, but now we need the analog of E for A1 instead of for A.

By definition, (dropping the 1’s for a moment)

E’A = I

and so

.

I am using Cohen’s

Ma = A’A

Me = E’E

and it turns out that Ma and Me are inverses. I will, of course, apply that equation for A and E to A1 and E1.

Here is Ma and its inverse Me:

Here is E1:

Now we do a Cholesky decomposition of Ma:

That is, the Cholesky decomposition has given us Ga such that Ga Ga’ = Ma.

Cohen had showed us that F2 = E Ga was an orthonormal matrix:

I want to point out now, that the lower-triangular structure of Ga says that the 3rd column of F2 is proportional to the third column of E1: in fact, I claim we have (although it’s difficult to make it clear)

(i.e. the 3rd column of F2 is .322916 times the 3rd column of E1).

Cohen also told us (I’ll confess that I haven’t figured out why) that the 1st column of F2 is proportional to the first column of A1: in fact, we have

Anyway, we have just applied his method to (f, E2, E3). And we have confirmed the realtionships between F2 and A1 and E1.

Fine, but what did we get? Did we get the same answers as from Gram-Schmidt?

Yes. The first column of F2 is equal to f1, the second column to f2, and the third column of F2 is equal to f3. (I could put up a drawing… but it looks exactly like the previous one.)

Lovely! And I’ll point out that we matched my f2 — not the negative of it. Even his method would require that we take the negative of his answer in order to get Figure 20.

The cross product

There is one other alternative computation I want to do.

We just saw that f1 was proportional to the first column of A1 — that’s how we chose it, a unit vector proportional to the fundamental f.

But we also saw that f3 was proportional to the third column of E1 — and so we could have computed f3 as soon as we had E1, just by normalizing the third column of E1. Cohen seems to have gone to an awful lot of work to get the matrix F2, given that both its first and third columns are known as soon as we compute E1 from A1. (This would work in all of his cases: for F2, normalize the first column of A and the third column of E; or — for F1 — normalize the first column of E and the third column of A.)

Anyway, that leaves f2, the second column of F2, to be determined.

Gee, it’s too bad we can’t compute a cross product for vectors of length 16.

Gee, who needs to?

Those vectors lie in a 3D subspace, so there is a change-of-basis transition matrix which leaves them with 3 nonzero components. In fact, the transition matrix is precisely v, and we know that new components y are related to old components x by

x = v y

so

since v is orthogonal. For old components f1, we get

t1 = v’ f1
= {-0.913557, -0.393456, -0.102984, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}

so we take y1 to be the three nonzero components…

y1 = {-0.913557, -0.393456, -0.102984}

For old components f3, we get

t3 = v’ f3
= {-0.361692, 0.901754, -0.236681, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}

and, again, we take y3 to be the three nonzero components…

y3 = {-0.361692, 0.901754, -0.236681}

Now, what we have are the i and k unit vectors (first and third columns), we want the j unit vector (second column). The right-hand-rule shows me — you can’t do this with your right hand on the mouse! — that

j = k x i,

so I compute

y2 = {-0.18599, 0.178973, 0.966114}

(which is a unit vector) and use

x = v y

to get the corresponding f. Oh, I need to pad y2 out with zeroes.

t4 = {-0.18599, 0.178973, 0.966114, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}

t2= v t4
= {0.013157, 0.141101, 0.24887, 0.137367, -0.126939, -0.342318, -0.483839, -0.440529, -0.221371, 0.106281, 0.346616, 0.341736, 0.187908, 0.0688953, 0.0186619, 0.00440436}

Quickly recall f2

f2 = {-0.013157, -0.141101, -0.24887, -0.137367, 0.126939, 0.342318, 0.483839, 0.440529, 0.221371, -0.106281, -0.346616, -0.34,1736, -0.187908, -0.0688953, -0.0186619, -0.00440436}

and we can confirm that t2 = – f2.

and that — presumably — is why he took the negative of f2: both his method and my Gram-Schmidt delivered left-handed coordinate systems; the cross product let us construct a right-handed coordinate system.

I think that’s pretty slick. Sure, it’s just a cross-product, but it was in a 3D subspace of .

Color: Cohen Figure 20, “an intriguing F matrix”

I struggled to make sense of pages 94-101 of Cohen’s “Visual Color and Color Mixture”.

I understand them now. Mostly.

The challenge was to derive a drawing of three basis vectors (Figure 20, p. 95), shown here :

He calls these three vectors both “an intriguing F matrix” and “This canonical orthonormal configuration….” But as far as I can tell, he never said where he got this particular one. He made it sound important, so I had to derive it. Figuring it out, as usual, was very educational.

What follows is a detective story.

I eventually decided that the most likely candidate was a “transformation matrix” on p. 101. There are five transformation matrices on that page, but the key is that exactly one of them is associated with an equal-energy illuminant — and that’s related to “EE” in the figure.

Let us proceed somewhat as I did, but without all the wrong turns and blind alleys. If you want a chance to wander those byways, then close this page, open a copy of Cohen, and go to work. OTOH, I’m about to take you down one, just one, wrong turn. Hey, it’s interesting and valuable.

First, let’s get that particular transition matrix of his. He denotes it but I will just refer to it as “TE” (outside of LaTeX).

He refers to the five T’s as “direct transformation matrices from the 1964 CIE X,Y,Z color-matching functions to orthonormal F matrices.” TE is specifically related to the equal-energy illumninant; the other four are related — in the same but as yet unexplained fashion — to the standard illuminants A, B, C, and D_65. (I have done the work to confirm that.)

OK. TE relates an A matrix to an F matrix. But how? “Transformation matrix” just doesn’t tell me enough.

A has 3 columns, so we could compute A TE, i.e. post-multiply A by TE. Since , , and are all the same size, we could use any of them in place of TE.

Well, A TE doesn’t look right — I dismissed it immediately when I first saw it — but eventually I saw the light. A TE does not produce an orthonormal matrix — but it’s only off by a factor of (and some change).

Let’s compute A TE and call it F….

Oh, first I had better get the A matrix.

Cohen uses the 1964 CIE tables for this, so I take the A matrix to be the 1964 XYZ (i.e. xyz bar) tables at 5 nm intervals. Not at 20 nm. My only target is a drawing, and I want a lot of detail. I plan to repeat part of these calculations, and to do some alternatives, at 20 nm in the next post. But for this one, I need detail.

I have found a set of these tables online. Let me show you, for checking purposes, a tiny subset of the XYZ tables as I store them:

My A matrix has 3 columns and 61 rows. The six rows of A corresponding to that part of the XYZ64 table are:

Mathematica aside.
There, I’ve finally shown you how I’m getting different A matrices whenever I change things. The table XYZ64 — like my other three tables — has wavelength for the first argument, and the Interpolation command lets me select any subsets, or even to interpolate. (It always gets the actual points right. In fact, I am usually simply getting subsets of the data, not truly interpolating. If I ever ask for 1 nm intervals, then I would be interpolating from the 5 nm table.)

For example, the following command gives me the entries for an A matrix from 400 to 500 in 20 nm intervals (and displays it nicely):

Table[Interpolation[XYZ64,i],{i,400,500,20}]//MatrixForm

So, once having set up my XYZ64 table, and the others, I can get any plausible A matrix I want.

end Mathematica aside.

For this problem, my A matrix is the entire 1964 CIE XYZ table at 5 nm intervals. Now we compute

F = A TE

and plot the results. Let me put the two drawings side-by-side (mine on the left):

That black curve could be his “EE” on a different scale, but my other two are utterly different from “I” and “II”.

Let me show you the first six rows of F:

Cohen said that F would be orthonormal. It is not. We would conclude that F was orthonormal if F’F = I. So we compute F’F:

Round it off:

I have no idea why his F — equivalently, his TE — is too large by . (The off-diagonal terms are almost zero, so these vectors are almost orthogonal. But their squared lengths are 2, so their lengths are . We could have fixed that by dividing TE by .)

Since the book was published posthumously, perhaps the editor was editing without computing. Perhaps he took what he found in Cohen’s notes or drafts and hoped that everything was consistent.

I assure you that the other possibilities (, , and ) do not come close to delivering an orthonormal F. Feel free to try them.

We conclude that there is a minor error in TE; it turns out that the other four transformation matrices on p. 100 also need to be divided by . Whatever he had in mind, he was consistent. Under the circumstances — posthumous publication, and the ease of computing unit vectors — this is rather unimportant — conceptually and computationally — once we know it.

More importantly, I think, we also conclude that the F produced by TE is not the F matrix in Figure 20. They seem to have — and actually do have — proportional first columns, but that’s all.

That is a little more troubling. I could have searched for a derivation of TE; instead, I searched for a derivation of Figure 20. I don’t care to do both. Let me be explicit: I don’t know how he got his TE matrix, nor any of the other 4 transformation matrices on p. 100. I do know how he got Figure 20.

It will be more than enough work to get Figure 20. I don’ t know where TE came from. I have some ideas, and I expect that you will too by the time we finish, but I’ m done with TE. Maybe I shouldn’t be, but enough is enough.

I should emphasize, however, that I am keeping the first column of F computed from TE, at least for a little while.

Oh, let me do one last thing with TE before I abandon it.

We know the relationship between A and TE: F = A TE. How would we express that in my usual terminology? Using T for TE, we have

F = A T

and we transpose, getting

F’ = T’ A’.

That says that each column of F’ is obtained by applying the matrix T’ to the corresponding column of A’. That tells us two or three things.

First and second, T’ is the transition matrix, and therefore his “transformation matrix” TE is actually an attitude matrix; third, he takes A as the new components and F as the old.

Just so I know.

Now let’s find a way to match Figure 20.

First, where did the first column — the black curve or “EE” — come from? (Yes, we showed that F produced from TE has a first column which is proportional to the “EE” curve. But where did this common vector come from?)

He tells us that it is proportional to the fundamental of the illuminant.

Ok. To get the fundamental, we need the projection operator. To get that, we need the SVD (Singular Value Decomposition) of the color-matching functions, i.e. of the matrix .

A’ = u w v’.

More to the point, we need the first three columns of the matrix v. (We need three columns because A’, like A, is of rank 3.)

Let me show you the first six (of 61) rows of v1:

Those 3 column vectors are an orthonormal basis for the preimage of the range of A’ (“the fundamental color space”). We form a projection operator by

R = v1 v1′

Now to get compute the fundamental of an equal-energy illuminant ee, we need to define the illuminant (I used all 100’s):

ee = {100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100}

Its fundamental f is the projection under R:

f = R ee =

{5.53525, 12.6614, 24.9595, 41.9659, 61.9674, 81.6212, 98.8882, 114.647, 125.821, 130.294, 129.298, 124.753, 116.918, 107.167, 95.2578, 80.7352, 68.1596, 60.0607, 54.8153, 53.7962, 55.5973, 59.7538, 65.4578, 72.9488, 81.1471, 88.5547, 95.6709, 102.8, 109.309, 114.255, 118.44, 122.816, 126.542, 129.052, 130.251, 129.621, 127.996, 126.529, 123.888, 119.146, 112.759, 105.059, 96.1651, 86.301, 76.0275, 65.9399, 55.8152, 45.5545, 36.3286, 28.6452, 22.185, 16.823, 12.5322, 9.19186, 6.64563, 4.72835, 3.33391, 2.32954, 1.61867, 1.12825, 0.779473}

Let me graph a scaled version of f along with the first column of F. (Scaled? I have to divide f by 470. Trust me.)

As I have done before, the black dots are one vector — in this case, the fundamental f — and the red curve is the other vector, in this case the first column of F.

We see that they are the same. Where we used to know that his TE gave us a column of F proportional to the EE curve in Figure 20, now we know what the curve is. Now I can dispense with F computed from TE, because I can use the fundamental f of the equal-energy illuminant instead.

The horizontal axis is off — I didn’t bother showing the wavelength.

That scaling by 470? Ha! The 1964 XYZ table at 1 nm intervals is over the range 360 to 830 nm. And 830-360 = 470.

(Probably “the length of the interval of integration”. Not that an integral is appropriate, in my opinion. And Cohen has pointedly been doing linear algebra, not calculus. But the literature is littered with integrals, even though the data is discrete. This is the mathematics I am skeptical of in all such works: they use calculus for their theory but linear algebra for their computations. (The example from Wyszecki & Stiles was closer to linear algebra, and the example they did two pages later is clearly not done by integration — no, I didn’t work that example.) I suspect it’s because a long time ago we pretty much restricted ourselves to continuous systems…. But I keep coming back to the fact that their data — the color matching functions — are discrete, and linear algebra is the appropriate tool. That’s my opinion.)

That’s where the 470 seems to be from, but I don’t really care. It would have made more sense for him to have used 400-700 instead of 360-830, and that’s one reason I don’t care about the scaling. Further, and more importantly, we will ultimately make a unit vector out of the fundamental f. Its scaling along the way just isn’t all that important.

And I didn’t scale the first column of F, so the in F is consistent with the 470 in f.

Sheesh! And I thought Fourier transforms had too many different scalings!

What interests me is that the first column of F is proprtional to f. And both appear to be proportional to the curve “EE” in figure 20. This is progress.

Now let’s get the dual basis E’. (I think I like idea of calling the rows of E’ a dual basis for the columns of A… and the columns of E a reciprocal basis for the columns of A… and, finally, the columns of E are a dual basis for the rows of A’. Row vectors and column vectors live in vector spaces which are dual to each other.)

Oh, say it compactly! The columns of E are simultaneously a reciprocal basis for the columns of Acolumns of A and a dual basis for the rows of A’. Numerically the same, conceptually different.

I recall the required computation as

E’A = I,

,

where doesn’t exist, but led me to the appropriate pseudo-inverse.

Like everything else in this example, E is too big to display, but here are its first six rows:

And here is a graph of the three rows of E’ (equivalently, the three columns of E):

NOTE that these are the “CIE 1964 XYZ fundamentals” on p. 163 of Cohen. Some of those resemble Figure 20… wait, E is not orthonormal, so none of them should be curves in Figure 20 (which are orthogonal, although not unit, vectors!) The resemblance is just that: not much more than vague similarity.

OK, let’s try getting an orthonormal basis, somehow.

I don’t know about you, but I want something I understand. Let’s apply Gram-Schmidt to a non-orthonormal basis. (This is almost certainly not what Cohen did. I’ll show you what I think he did in the next post. We’ll get the same answer.)

The first two non-orthonormal bases I started with were

f, E1, E2

and

f, E1, E3.

In retrospect, my third choice

f, E2, E3

seems more obvious. OTOH, more obvious or not, what matters is that it does, in fact, generate Figure 20!

As is my custom, I am going to describe Gram-Schmidt for this 3D case. The general method just extends it from 3D to n-dimensions.

We have our first vector f. Normalize it, by

.

This gives us our first unit vector, f1, for the new basis.

f1 = {0.00827079, 0.0189187, 0.0372946, 0.0627057, 0.092592, 0.121959, 0.147759, 0.171307, 0.188002, 0.194686, 0.193198, 0.186406, 0.174699, 0.16013, 0.142335, 0.120635, 0.101844, 0.0897429, 0.0819053, 0.0803825, 0.0830737, 0.0892844, 0.0978074, 0.109, 0.12125, 0.132319, 0.142952, 0.153603, 0.163331, 0.170721, 0.176973, 0.183512, 0.189079, 0.192829, 0.194622, 0.19368, 0.191252, 0.18906, 0.185114, 0.178029, 0.168484, 0.15698, 0.14369, 0.128951, 0.113601, 0.0985278, 0.0833993, 0.0680677, 0.0542824, 0.0428017, 0.033149, 0.025137, 0.0187257, 0.0137345, 0.00992992, 0.00706512, 0.00498154, 0.00348081, 0.00241862, 0.00168584, 0.00116469}

Normalize E2 and E3, as we did f, just to start with vectors of the same size. Call the unit vectors e2 and e3.

e2 = {-0.00326355, -0.00739089, -0.0142972, -0.0235132, -0.0336395, -0.0422167, -0.0480022, -0.0514717, -0.0506656, -0.0450514, -0.0346783, -0.0212322, -0.00302856, 0.0178248, 0.0440051, 0.0714143, 0.0954016, 0.120945, 0.141852, 0.165312, 0.189969, 0.213125, 0.235048, 0.255043, 0.270026, 0.276282, 0.276115, 0.273313, 0.264794, 0.24833, 0.226531, 0.201183, 0.171297, 0.136324, 0.0977252, 0.0577811, 0.0185254, -0.0185732, -0.0524476, -0.0805823, -0.102238, -0.116778, -0.124173, -0.124757, -0.11956, -0.109979, -0.0974696, -0.0829623, -0.0686652, -0.0558823, -0.0443658, -0.0340783, -0.0255961, -0.0189148, -0.013768, -0.00982428, -0.00696574, -0.0048783, -0.00340378, -0.00233611, -0.0016479}

e3 = {0.0134212, 0.0307963, 0.0609534, 0.103042, 0.152976, 0.202366, 0.245956, 0.285754, 0.313861, 0.325009, 0.321797, 0.308947, 0.28685, 0.259425, 0.22519, 0.182631, 0.143994, 0.114635, 0.0919794, 0.0761547, 0.0649611, 0.0570429, 0.0503958, 0.0457693, 0.041089, 0.0352787, 0.0289087, 0.0223688, 0.0152318, 0.00702492, -0.00182114, -0.0115187, -0.0217216, -0.0324157, -0.0433804, -0.0537281, -0.0633608, -0.0725008, -0.0801518, -0.0851463, -0.0875296, -0.0873093, -0.08454, -0.0794315, -0.0725582, -0.0646114, -0.0558614, -0.0465045, -0.0377562, -0.0302361, -0.0237077, -0.0180942, -0.0135352, -0.00996527, -0.00722962, -0.00515145, -0.0036426, -0.00254819, -0.00177437, -0.00122704, -0.000856806}

Now apply Gram-Schmidt to the pair f1 and e2. The vector e2 is not orthogonal to f1, but it can be written as the sum of something parallel to f1 plus something perpendicular to f1. We will find the parallel part using the dot product, and get the perpendicular part by subtraction.

The component of e2 wrt f1 is… (since f1 is a unit vector), so we make a vector of that length, parallel to f1, and subtract it from e2. The result of that vector subtraction is perpendicular to f1. We normalize it and call it f2.

.

Then

f1, f2

is an orthonormal basis that spans the same subspace as

f1, e2.

f2 = {-0.00662133, -0.0150664, -0.029409, -0.0488852, -0.0710293, -0.0913202, -0.107276, -0.119901, -0.125361, -0.121889, -0.110234, -0.0932842, -0.0693866, -0.041572, -0.00683967, 0.0306849, 0.063449, 0.0953437, 0.120668, 0.14633, 0.171675, 0.194084, 0.214299, 0.231442, 0.242825, 0.245323, 0.241118, 0.234089, 0.221296, 0.200893, 0.175217, 0.145636, 0.111572, 0.072756, 0.0308042, -0.0115508, -0.052607, -0.0914459, -0.126173, -0.153574, -0.173117, -0.184308, -0.187183, -0.182227, -0.170858, -0.154906, -0.135801, -0.114483, -0.0939761, -0.0759605, -0.0599912, -0.0459572, -0.0344598, -0.0254257, -0.0184818, -0.0131801, -0.00933457, -0.00653426, -0.00455539, -0.00313629, -0.00220307}

Now what do we do with e3? We want to break off the piece which lies in the plane spanned by f1 and f2. That’s just the vector

.

When we subtract that from e3, what remains is perpendicular to the plane spanned by f1 and f2; we will normalize it to length 1 and then call it f3. (It’s crucial that f1 and f2 are an orthonormal basis for that plane; otherwise we would need an explicit reciprocal basis to find components of e3 wrt f1 and f2.)

f3 = {0.0108058, 0.0248259, 0.0492157, 0.0833774, 0.124048, 0.164396, 0.200091, 0.232718, 0.255782, 0.264988, 0.262322, 0.251597, 0.233086, 0.2101, 0.18122, 0.145023, 0.111701, 0.085438, 0.064457, 0.0480153, 0.034675, 0.0232252, 0.0120345, 0.00174899, -0.00919785, -0.0209678, -0.0332539, -0.0457821, -0.0586168, -0.071629, -0.0848832, -0.0992835, -0.113845, -0.128148, -0.141854, -0.153559, -0.163729, -0.173431, -0.180567, -0.183118, -0.181455, -0.175821, -0.166361, -0.153478, -0.138238, -0.121869, -0.10453, -0.0863848, -0.0696761, -0.0554858, -0.0433135, -0.0329815, -0.0246352, -0.0181132, -0.0131248, -0.00934719, -0.00660276, -0.00461709, -0.00321259, -0.00222783, -0.00154979}

(Maybe I shouldn’t have printed all those, but if you have the 1964 xyz bar tables and you’re trying to verify my calculations….)

So, it’s time to look at what we got. Here are f1, f2, f3. They are an orthonormal basis, and they span the same space as E1, E2, E3 — the fundamental color space, the pre-image of the range of A’.

The brown looks like Cohen’s “I”. But the mauve curve is off. Still….

Do I need the negative of the mauve curve?

Yes!

And multiply by 44 ! (No, I have no idea why! I just kept changing the multiplier until I matched his scale on Figure 20.)

Having found them, make them pretty — i.e. set the horizontal scale.

Here, let me put them side-by-side. They look very close, but not exact. Nevertheless, I am content: his numbers tend to be off a little, so what’s the point of trying to match them perfectly?

So. Start with the non-orthonormal basis f, E2, E3 and apply Gram-Schmidt to get the orthonormal basis f1, f2, f3. Then take the negative of f2.

I must emphasize that his scale (i.e. multiply by 44) does not represent orthonormal vectors. Pairwise orthogonal, yes — but the f1, f2, and f3 I derived were specifically unit vectors. The scale on his figure requires that I multiply by 44 (or something very close to it). Once we have unit vectors, we can scale them any way that is convenient. Here it was convenient to match his — frankly, incomprehensible — vertical scale.

Oh, let’s close with one last calculation, to compute the “transformation matrix” for Figure 20. We have found the F matrix for Figure 20, and we have the A matrix. What is the transition matrix T between them? (Actually, we know we want the attitude matrix instead.) Recall that he takes A as the new, and F as the old, components.

So, as usual, we compute each column of F’ (old) by applying T to each column of A’ (new). That is,

F’ = T A’

F = A T’

and we write the false but suggestive

which gives us the true

using the appropriate pseudo-inverse.

So we compute and display T’, and remember that, as the transpose of a transition matrix, it is an attitude matrix…

We can compare that to the transformation matrix TE we tried at the beginning:

As we would expect — since the F produced by TE did not match Figure 20 — these are different.

Real Short Summary

I matched his Figure 20 by applying the Gram-Schmidt orthogonalization procedure to the vectors f, E2, E3… where f was the fundamental of an equal-energy illuminant, and E2 and E3 were two of the basis vectors dual to the rows of A’.

Summary

We really started with two separate items: a “transformation matrix” TE — one of five on p. 100 — and Figure 20 on p. 95 of Cohen. The transformation matrix TE turned out to be an attitude matrix describing a new A matrix in terms of an old F matrix. (Yes, that’s backwards, but so long as we knew it, it didn’t matter. This is why I carefully distinguish attitude and transition matrices; so long as I know what I’ve been handed, I can use it properly and confidently.)

We saw that F = A TE reproduced only the “EE” curve of his Figure 20. And his F matrix failed to be orthonormal by a factor of .

We then saw that the “EE” curve was (proportional to) the fundamental f of an equal-energy illuminant. In order to find the fundamental of a spectrum, we had to find a projection operator R onto the pre-image of the range of A^T. The fundamental f is then the result of applying the projection operator R to a constant vector (“equal energy”).

By trial and error, I found that the F matrix for Figure 20 could be (very nearly) reproduced by applying the Gram-Schmidt orthogonalization procedure to the vectors f, E2, E3, where E1, E2, E3 are a basis dual to the rows of A’. We found the dual basis as the rows of and then transposed E’ to get E, a matrix of 3 columns.

In the next post, I will repeat these calculations quickly using 20 nm intervals, and I will show you a few alternatives. In particular, I will show you that his method (a Cholesky decomposition) gives the same answer — even to having to switch the sign of the second column!

Transition matrix: to be or not to be

Cohen (“Visual Color & Color Mixture”, see the bibliography) did something very interesting. In fact, he did something useful which I had never seen before.

Although this post uses some matrices which we saw in the color posts, I think this can stand on its own: you need not have read the color posts. But if you are specifically interested in color, or in Cohen’s work, this post is very relevant.

He was trying to describe how to find a transition matrix between two given data matrices. This will come in handy — very handy! — whenever people give the alleged result of an unspecified linear transformation of a data matrix.

On p. 93, he wrote

D T = C in general

F T = A in particular

with

.

Because D’D (that is, ) is invertible, we infer that D and C are taller than wide, like the A matrix… like our usual data or design matrix X.

We also infer that T is a small matrix: its size is cxc, where c is the number of columns of D. D and C are the same size.

I do have to point out that in his equation (and his notation)

D T = C,

T is an attitude matrix! I will change that. (Trust me: you want me to change that.)

Since we have both XYZ (i.e. 1931 xyz bar) and RGB (1931 rgb bar) tables, and since they are supposed to be related by the transformation matrix T31…

… we should try the idea embodied in that equation for T. Assuming that there is a transition matrix T, let J and K stand for the XYZ and RGB matrices respectively.

As usual, let me use 20 nm intervals instead of the published 5 nm intervals. We have

In case you have not read the color posts, all I am saying is that J and K and T31 are given matrices, and T31 is supposed to be a transition matrix between J and K.

We do know — or take it as given — that for the transition matrix T31, old is XYZ ~ K, so let me write

K’ = T J’

That says T is a transition matrix which, applied to each column vector of J’, delivers the corresponding column of K’.

Transpose:

K = J T’

and then pretend exists (even though J is rectangular)

and then substitute the appropriate pseudo-inverse:

.

(The inappropriate pseudo-inverse would use JJ’ instead — which cannot be of full rank, given the assumed shape of J, hence cannot be invertible.)

That equation defined the transpose T’, but here’s T:

Recall the given T31:

By doing it my way, I constructed the counterpart to T31, which is the published transition matrix between XYZ and RGB (1931) with XYZ considered the “old” basis.

T and T31 are very close, but not the same.

Even if I used the full tables, at 5 nm intervals, we would not get T31.

First, I am certain that if there were a transition matrix T, then we must have computed it. Since this computation did not yield T31, I infer that T31 is not the transition matrix.

Second, is T the transition matrix? We should apply T to the appropriate table and see what we get. We had assumed

K’ = T J’, or equivalently

K = J T’

so let’s compute a counterpart to K, namely Q = J T’:

and compare it to K, via Q – K:

(I keep using pictures for things like that, namely scientific notation, just because it’s easier than editing the LaTeX.)

The largest entry in absolute value is… 0.000069765 . The three columns of Q look like this:

Pretty small, and pretty random, but still not really zero. Well, maybe it wouldn’t be. Let’s try our magic formula on J and Q, since they are related by the transition matrix T. From

,

I compute what I call t’:

,

and then display the transpose of the transpose, t” = t:

Recall T:

So. That equation for T’ does recover the transition matrix if there is one. The discrepancy between T and T31 is real.

Maybe I’m making too big a deal of this. I infer that the XYZ and RGB tables do not span exactly the same subspace. They almost do, but not quite.

To be more explicit: the 1931 XYZ (xbar, ybar, zbar) and RGB (rbar, gbar, bbar) tables are not exactly related by any transition matrix whatsoever, although the published T31 is close to being one for them.

I have found a marvelous example — disappointing as a published result but marvelous as a bad example — where the discrepancy is significant. I expect to show it in another post, soon.

We’ve been looking at the computation in one way: if there is a transition matrix T such that

K’ = T J’

then

.

What if there is no such transition matrix T? We can still define T by that equation, and compute (still assuming J is of full rank, with more rows than columns), but what are we getting, if T itself does not exist as a transition matrix?

That is, what if we compute T, and find that

?

Well, this will be clearer in the next post, but…. That equation for T’ looks an awful lot like the normal equations for a least-squares fit.

What I will illustrate in the next post is that generates a transition matrix to the subspace spanned by J. That is, T is a transition matrix — any invertible matrix is! — but it relates J to a subspace other than (the one spanned by) K.

(In our case, we had 21 observations; the 3 columns of J and of K define 3D subspaces of . Those two 3D subspaces are almost the same, but not exactly.)

Let me try to explain here it without a good (I mean, bad) example. If J and K span the same subspace, then any one column of K (think of it as the dependent variable y) can be written as a linear combination of the columns of J (think of it as the X matrix).

But if J and K do not span the same space, then the best we can do is to find that linear combination (call it yhat !) of columns of X which is as close to y as possible.

They really are examples of y and yhat and X.

And what is T’ ? Well, if we write it with y and X…

we might recognize that as in ordinary least squares.

Let me show you that, while providing another example of the utility of the pseudo-inverse.

The derivation of the normal equations for ordinary least squares is fairly complicated, if only because it involves the derivative of a matrix. But the pseudo-inverse would let us recover the equations themselves.

Deriving them is one thing; recovering the form of them once we know they exist is another.

We write our model (this is a quick & dirty recollection of the answer, not a derivation):

.

(I guess I have to talk about that. True is . But I write y instead of yhat, as though we have equality rather than a projection onto a subspace.)

We pretend that X is invertible, and write

.

Since X is not invertible, but X’X is (by assumption X is of full rank with more rows than columns), we replace the non-existent inverse by the pseudo-inverse and write

.

The real derivation gave us that very answer. Once I know that, I never need to compare them again. I have a plausibility argument that gives me the known right answer.

This is analagous to confirming in the previous color post, among others, that if I compute

then I have a dual (or reciprocal or biorthogonal) basis E such that

E’A = I.

In that case, I prefer to verify it computationally every time.

Incidentally, this trickery with the pseudo-inverse would also generate a true result if I had started with yhat instead of y — it is true that , but it’s not very useful because we don’t know yhat!

You might have noticed that corresponds, algebraically, to a column of T’. Asking if the first column of K is a linear combination of all the columns of J gives us the first column of T’; asking if the ith column of K is a linear combination of all the columns of J gives us the ith column of T’. This, too, will be in the next related post.

Summary

Given two data matrices J and K, tall and thin, we can define and compute

.

If the individual observations are in fact related by a transition matrix M, so that

K’ = M J’,

then our computed T is M:

T = M.

This is pretty slick. Any time someone hands us some data and the alleged result of applying an unspecified linear transformation to it — we can find that linear transformation if it exists, or demonstrate that it does not exist (in which case, they goofed).

And if there is no such matrix M? Well, it’s simple enough: if there is no such M — and J and K are both of full rank — then our computed T must be something else, so we do not have equality:

.

That’s all it takes for us to know that there is no transition matrix M between J’ and K’. If there were, it would be T, but that didn’t work.

(That’s if J and K are both of full rank; if only one is not, then one of J’J or K’K is actually not invertible — as JJ’ and KK’ are not invertible — and T can only be computed in one direction. I expect that T is not invertible if J and K are not both of full rank.)

And when there is no transition matrix M, we nevertheless get a transition matrix T, but it’s between J’ and Q’ = T J’.

I will freely confess that I did not know any part of this before I read Cohen: neither that we could find the transition matrix so easily; nor what we were finding if it were not a transition matrix.

(Cohen implicitly assumes that there exists a transition M; he does not discuss what happens if there isn’t one.)

Next post on this topic, a numerical illustration.

Color: decomposing the RGB A transpose matrix

(abuse of) terminology

Sometimes I get tired of writing “xbar, ybar, zbar tables” — and I just write “xyz bar tables” or even “XYZ tables”. Similarly for rbar, gbar, bbar tables — rgb bar. I’m not talking about anything new, just abbreviating the names.

Introduction

This is the third post about the example on p. 160 of W&S. Once again, I am going to decompose the reflected spectrum into its fundamental and its residual.

This time, however, I’m going to use the rbar, gbar, bbar tables (RGB) instead of the xbar, ybar, zbar (XYZ) tables. They did not do this.

I’ll tell you now there is one little twist in these calculations. We will need the ybar table, because we still need to use it to scale our results.

In addition to showing, as I did previously, the dual basis spectra and the orthonormal basis spectra for the non-nullspace, I will display the orthonormal basis for the nullspace.

As usual, I am using the CIE 1931 tables, and I am working at 20 nm intervals.

Our ingredients were: a reflectance spectrum, an illiminant spectrum, and some choice of color matching functions. For this post, I am choosing a different set of color matching functions.

Our results in the previous post were a few sets of 3 numbers. The final result was called xyz — meaning the triple (x, y, z) — for plotting a point on the CIE chromaticity chart. We will not do this again.

That result came from a triple denoted XYZ, which are values in XYZ color-space. That space has the advantage of being device-independent. We could — and did — convert those values to an RGB color space, specifically the one defined by the rgb bar (i.e. rbar, gbar, bbar) color-matching tables.

This post will compute the RGB values directly, and we will find that we get almost exactly the same RGB values. (We really ought to.)

The XYZ result in turn came from a triple for which I have no convenient notation, but I keep saying without proof, and confirming by computation, that that triple is the components of the fundamental of the reflected spectrum wrt the basis dual to the chosen color-matching functions.

Let us go see that, this time for the rbar, gbar, bbar tables instead of xbar, ybar, zbar.

The same old calculations, but with a twist

But first, let me review the calculations. Feel free to skip to the next section.

Our first computation is to get the reflected spectrum, which is the pointwise product of the illuminant spectrum and the reflectance spectrum. For this problem, the reflected spectrum is the perceived spectrum, the input to the eye of “the standard observer”. The data itself is in the fourth case of the first post about this example.

Here is what the three of those spectra looked like:

As before, I’ve decided to avoid red, green, and blue — except where they are appropriate. For this post, they will be appropriate in some cases. The illuminant is shown in white, the reflectance spectrum is shown in purple, and the reflected spectrum is shown in yellow.

Having said that, I am going to use red, green, blue colors for the three columns of our A matrix, because these entires are rbar, gbar, bbar values. These are the 1931 color matching functions (i.e. a subset, at 20 nm intervals)

BTW, we can tell that we are using RGB (i.e. rbar, gbar, bbar tables) because one of these functions is clearly negative over part of its range. In fact all of them are negative somewhere, but it’s only obvious for the red one.

Since I’ve just introduced new data, let me display it:

Our second computation (or third, this order doesn’t matter) should be to apply the matrix A^T — the (transpose of) the color matching functions A — to the reflected spectrum S:

.

As before these numbers are components of (the fundamental of) the reflected spectrum wrt the basis dual to . In Cohen’s notation, these 3 numbers are the components wrt the columns of the E matrix.

The other computation, whether we do it second or third, is to find the constant for normalizing things, using the illuminant spectrum; that is, to set the white point (1,1,1) in our color space.

Here’s the twist. We still have to use the second column of the XYZ (xyz bar) tables; not one of the columns of our new A is the total cone response. (Nor is the sum of them, and I’ll confess that I don’t know why. My WAG (my guess) is that the three color matching functions each have to be weighted differently in order to sum to the total cone response.)

The number we got by applying the second column of the XYZ tables to the illuminant was 529.764 .

Now we move on to the fourth computation: we divide our 3 values by that to get RGB.

RGB = {0.133246, 0.115616, 0.0903067}.

It is encouraging to see that these RGB values are in [0, 1]. These are tristimulus values in an RGB color space, namely the one I keep showing on the chromaticity chart.

Let me be clear. The red, green, and blue dots on this diagram are the projections of (1,0,0), (0,1,0) and (0,0,1) onto xy-space. The values I just got are in this RGB space.

Oh, in the previous post we transformed the computed XYZ values into RGB; what did we get?

{0.133246, 0.115616, 0.0903103}

So. Almost exactly the same answers: only the B value differs. That’s good. I said I would confirm this equivalence.

The blue value is off a tiny bit. (The simple fact is that the two sets of tables — xyz bar and rgb bar — and the transition matrix between them are not consistent. In fact, it appears that the columns of the two tables do not quite span the same 3D space. If I’m right, I could fix this by using one table and the transition matrix to construct a very slightly different version of the second table. But bear in mind that I am printing far more digits than W&S showed.)

the dual basis

Now, let’s get the dual basis E. We recall that instead of

,

we use the pseudo inverse and write

It is easier to work with but easier to display E:

Check that we do have . I didn’t say it last time, so let me emphasize: I didn’t prove that the pseudo-inverse would give me the dual basis. This easy “check” is the real confirmation that my recipe actually works. This check is crucial, not merely convenient.

Let me call the basis vectors E1, E2, E3, and let’s look at them. In fact, let me juxtapose the dual basis (right) with the color matching functions (left).

If we multiply each of those dual basis vectors (E1, E2, E3) by our original three numbers…

{70.5889, 61.249, 47.8412}

and add the results.. we get numbers…

{0.243727,2.51476,23.848,64.4959,61.7521,33.9471,25.2246,
43.8826,66.214,82.1295,90.1689,83.8434,60.5088,30.3679,
10.993,3.10116,0.751257,0.192395,0.0458083,0.010994,0.}

and a picture…

where I have shown the sum in white.

Let’s look quickly at the reflected spectrum (in yellow) and that sum (still in white):

Yes, we’ve seen that before. And as before, that white curve is the fundamental. It differs from the yellow curve by a metameric black, a spectrum in the null space of .

As before, we could stop here. We have gotten RGB values from the reflected spectrum… we have found the dual basis… we have found the fundamental of the reflected spectrum… we know that the first three numbers we computed () are the components of the fundamental wrt the dual basis.

The SVD, orthonormal basis, and projection operator

Now let’s do the SVD (Singular Value Decomposition), because we can get a nice orthonormal basis for the entire domain of including its null space, i.e. for spectra including the metameric blacks.

As before, we will get an orthonormal basis v1 for the pre-image of , then a projection operator R onto that pre-image, then recompute the fundamental by applying the projection operator to the reflected spectrum.

First, let

.

(The first time I did this, I let . Now we’re more experienced, but if you want to use B instead of , feel free.)

As always, our vectors have length 21 because we used 20 nm intervals from 380 to 780 nm. Then, the leftmost 3 columns of v are an orthonormal basis for the pre-image of , and the rightmost 18 (i.e. 21-3) are an orthonormal basis for the nullspace of .

Here are the leftmost 3 columns of v, denoted v1 as is my habit:

OK, I have to know: did v1 change? Yes. I’ll show you, soon.

Let me juxtapose the orthonormal basis v1 (right) with the dual basis E vectors. Both bases span the same space, the pre-image of , but they are different.

And in this case, I can’t see that any vectors in one basis look similar to vectors in the other basis.

Let me recall an image from the previous post: here are the dual basis E and the o/n basis v1 from the previous post. This is the old counterpart to the preceding drawing.

And the right hand pieces of those two drawings are the new v1 and the old v1. We can see that they are quite different.

Having an orthonormal basis for the preimage of , we know that we can construct a projection operator R from the domain of onto the preimage; we just take the product of v1 and in the other order:

.

It’s 21×21, as before, so I’m not going to show it to you. But, for checking, I’ll show you the first column:

And I’ll tell you that the last column, as before, is identically 0.

the fundamental and the residual

So. Apply the projection operator R to the reflected spectrum S to get the fundamental f…

f = {0.243727,2.51476,23.848,64.4959,61.7521,33.9471,25.2246,
43.8826,66.214,82.1295,90.1689,83.8434,60.5088,30.3679,
10.993,3.10116,0.751257,0.192395,0.0458083,0.010994,0.}

and subtract to get the residual n = S – f…

n = {4.85627,26.2996,19.116,-14.6684,-7.32845,18.6715,27.0686,
15.2246,3.00319,-13.0295,-23.7795,-12.7434,13.5977,41.0282,
59.1236,63.6887,60.2519,52.106,59.2832,33.0258,39.625}

Did we get the same fundamental? They are very close. Again, I think the real problem is that the two tables are not exactly consistent. Here’s

old f minus new f =

{-0.00362232,0.000221807,0.00281357,-0.00101889,-0.00138667,
-0.00346164,-0.00173991,-0.00126245,-0.00325203,-0.000486578,
0.00271479,0.00422385,0.00120929,-0.00248099,-0.00107542,
0.0013427,0.00223839,-0.00306751,0.00309787,0.00376061,0.}

We already know that the fundamental is the same whether we compute it using the dual basis or using the projection operator. Let me illustrate that the fundamental computed using RGB is the same as the fundamental computed using XYZ.

I have plotted points using the previous fundamental, and I have plotted a piecewise linear curve using the new fundamental:

I should probably emphasize that although this drawing looks like one from the previous post, it illustrates something else. The apparently identical drawing in the previous post showed the equivalence of the fundamental computed two ways using XYZ; this drawing shows the equivalence of the fundamental from the previous post with the fundamental from this post.

And the null? Close, too. We have

old n minus new n =

{0.00362232,-0.000221807,-0.00281357,0.00101889,0.00138667,
0.00346164,0.00173991,0.00126245,0.00325203,0.000486578,
-0.00271479,-0.00422385,-0.00120929,0.00248099,0.00107542,
-0.0013427,-0.00223839,0.00306751,-0.00309787,-0.00376061,0.}

We can check that n is in the nullspace of by computing :

{-6.77791*10^-14, 6.72509*10^-15, -3.50254*10^-14}.

We should also check that the fundamental has the same 3 components as the entire reflected spectrum. That is, we had

.

for applied to the reflected spectrum. If we apply to the fundamental, we get

.

The same. As they should be.

You know what? There’s another pair of calculations I can do, to make a point.

What are the components of the fundamental wrt the orthonormal basis v1? We just dot f into the v1 matrix, because its columns are orthonormal basis vectors:

{-129.491,141.916,-79.7421,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}.

That confirms that no part of the fundamental lies in the nullspace. That was the whole point of the projection operator.

And what are the components of the residual? I hope you’ve got it: exactly the first three are zero:

{0,0,0,-24.6623,-15.3786,17.8955,33.3365,29.6755,19.4414,
-0.412828,-18.9652,-15.4907,8.59491,37.8598,57.8821,63.3301,
60.1642,52.0835,59.2778,33.0245,39.625}.

The decomposition looks the same as in the previous post. The reflected spectrum is in yellow, the fundamental in white, and their difference n, the residual, is in black.

NOTE that although reflected and fundamental are non-negative, the residual is not. The good news is that, in principle, we could generate the fundamental spectrum, as well as the reflected. But I can’t imagine how to generate a spectrum with negative values.

We face this problem not just for the residual, but for just about every basis vector we’ve computed!

(The middle vector of the orthonormal basis, i.e. the middle column of v1, is the only exception so far; it is everywhere non-negative. But that’s only for the orthonormal basis for the RGB color matching functions; for the previous post, every column of v1 has some negative values.)

Because the basis vectors are, frankly, weird spectra, don’t expect anything nice if you try to find their RGB, XYZ, or xyz values. Most of them lie outside the CIE chromaticity chart! Incidentally, I have computed Cohen’s F1 and F2 bases for this case, and some of them are still negative, too. Well, they have to be: things have to cancel out so that the dot products among the orthonormal vectors are actually zero.

But while we’re looking at the linear algebra, let’s take a look at the orthonormal basis for the null space.

the null space basis

Now for something new.

We have 18 basis vectors (metameric black spectra) for the nullspace. It’s rather confusing to look at all of them… it’s not much better to look at a few of them close together. In the top left, we have the 1st 5th, 9th, and 13th (1,5,9,13)… the top right is 2,6,10,14… bottom left is 3,7,11,15… bottom right is 4,8,12,16:

There are two more to show (17 & 18) and I thought I show a few more adjacent pairs, too. I have added 1 & 2 and 9 & 10:

While I’m playing around…. What about four spectra centered around 500? And add them up, to get another metameric black.

Let me do one last thing.

Graphs at 5 nm

What would the bases look like if I used the full 5 nm tables? The fundamental depended on the reflected spctrum, but the bases depend only on my choice of color matching functions. So, let me show them for the 1931 CIE rgb bar tables at 5nm intervals.

Here are the orthonormal basis vectors “v1″ — at 20 nm on the left, 5 nm on the right:

And the dual basis vectors “E”:

And the dual basis “E” juxtaposed with the orthonormal basis “v1″:

I might as well do this for the XYZ color matching functions, too. The more detailed curves are from the 1931 xyz bar tables, at 5 nm intervals.

Here are the orthonormal basis vectors “v1″:

The dual basis vectors “E”:

And the dual basis “E” juxtaposed with the orthonormal basis “v1″:

Color: decomposing the A transpose matrix and the reflected spectrum

This is the second post about the example on p. 160 of W&S. Here’s the first. (Incidentally, I am going to decompose the reflected spectrum into its fundamental and its residual. I do not think that they did this. Oh, this is what I did for Cohen’s toy example, but now it’s for real. Approximate, but for real.)

review the XYZ for this

Let’s review one case from the previous post: 1931 CIE at 20 nm intervals. As is my custom, I took a subset of their values — I did not average the numbers to account for the data I omitted. If this approximation isn’t good enough, simply go to 5 nm or 1 nm intervals.

Our ingredients were: a reflectance spectrum, an illiminant spectrum, and some choice of color matching functions.

Our results were a few sets of 3 numbers. The final result was called xyz — meaning the triple (x, y, z) — for plotting a point on the CIE chromaticity chart.

That result came from a triple denoted XYZ, which are values in XYZ color-space. That space has the advantage of being device-independent. We could — and will shortly — convert those values to an RGB color space, specifically the one defined by the CIE 1931 rgb bar (i.e. rbar, gbar, bbar) color-matching tables.

That XYZ result in turn came from a triple for which I have no convenient notation, but I keep saying without demonstration that that triple is the components of the fundamental of the reflected spectrum wrt the basis dual to the chosen color-matching functions. (Try saying that with a mouthfull of rocks!)

Let us go see that.

But first, let me review the calculations. Feel free to skip to the next section.

Our first computation is to get the reflected spectrum, which is the pointwise product of the illuminant spectrum and the reflectance spectrum. For this problem, the reflected spectrum is the perceived spectrum, the input to the eye of “the standard observer”.

Here is what the three of those spectra looked like:

For this post, I’ve decided to avoid red, green, and blue, except where it is really appropriate. The illuminant is shown in white, the reflectance spectrum is shown in purple, and the reflected spectrum is shown in yellow.

For this example, I had chosen the 1931 xyzbar tables for our color functions, and they looked like this:

Once again, I’ve changed the colors. X is orange, Y is white, and Z is yellow. It might help if you remember that the Y column of this A is nothing like green, being instead the “photopic response” — the total energy received by the cones — of the standard observer; white seems a sensible color for it.

BTW, we can tell that we are using XYZ (i.e. xbar, ybar, zbar tables) because these functions are never negative; we recall that each of the rbar, gbar, bbar color matching functions has some negative values.

Our second computation (or third, it doesn’t matter) should be to apply the matrix — the (transpose of) the color matching functions A — to the reflected spectrum S, getting:

I claim that these numbers are components of the fundamental of the reflected spectrum wrt the basis dual to . In Cohen’s notation, these 3 numbers are the components wrt the columns of the E matrix.

Please be patient. We’ll get there real soon now.

The other computation, whether we do it second or third, is to normalize things using the illuminant spectrum; that is, to set the white point (1,1,1) in our XYZ space.

We did this easily because the second column of A — the second row of — is the total cone response — the photopic response — of the standard observer.

Here’s the second row of applied to D65: 529.764

Now we move on to the fourth computation: we divide our 3 values by that to get X,Y,Z.

XYZ = {0.673539, 0.669432, 0.511754}

NOTE that the XYZ values are in [0, 1]. These are tristimulus values in a sort-of color space. (Let’s see. It is not HSI, even though that middle number is an intensity. In fact, it does not directly specify a color; that’s done by a change-of-basis. Maybe I should call it a pre-color space? Maybe I’m getting carried away here.)

And then our fifth computation is to get x,y,z. Here’s the sum X+Y+Z: 1.85472

and we divide XYZ by it to get xyz values: xyz = {0.363148, 0.360933, 0.275919}.

Where was that on the chromaticity chart? The yellow dot.

All that was done in the previous post. Now let’s move on.

Moving on.

Get the RGB

Let’s do one last thing before putting this down. I want to get RGB values from these XYZ. I have in my hands, and have shown you before, the transition matrix T31 (having once needed the 1964 CIE, I have renamed some things)…

It treats XYZ as old, RGB as new, so we have

XYZ = T31 RGB,

so

.

We’ll come back to these in the next post. These should be RGB color-space values using the 3 sources shown on the CIE chromaticity diagram; alternatively, these should be the RGB values computed directly from the rbar, gbar, bbar tables. That’s what I’ll confirm in the next color post.

Moving on, finally, there are two things we need to do. Actually, we might get away with only one of them, and that’s the dual basis, so let’s go get it first.

getting the dual basis E.

As I have said before, given a basis with transition matrix A we can find a basis with transition matrix E dual to it by requiring that the dot products of the vectors satisfy

.

That is, by laying out the E vectors as rows and the A vectors as columns, we compute all of the possible dot products as .

A digression on dual bases

I’ll remind you that we are working in finite dimensional vector spaces, and the dual basis and reciprocal basis look exactly the same. (So if you need to review, look at my posts on reciprocal bases. See the tag cloud.) But I think it’s worthwhile to think of reciprocal bases as living in the same vector space, while dual bases live in a vector space V and its dual V*.

Oh, the dual space V* is the vector space of all linear functionals whose domain is V, i.e. all real-valued linear operators from V to R.

I’ll also remind you that although one usually starts with a basis in V and finds “the dual basis” in V*, we are going backwards. The rows of are a basis for V*, and we wish to find a matrix E whose columns are a basis for V, and more, a dual basis to the rows of .

What a tangled knot! In my head I was going backwards, but I computed it the other way! What I really computed was a basis dual to A — but that is the same as finding a basis E dual to . The two equations

and

are equivalent. I chose to solve the first instead of the second, even though I was thinking of the second. (Aside: it’s a very good thing that the dual of the dual space is the original space, for finite dimensional vector spaces.)

Finally, I suppose I should say again that although we often compute the dot product of two vectors (u and v, say) by writing one as a row (say u) and the other (v) as a column, we may view that matrix product as the effect of applying the linear functional u to the vector v. And what I did was the equivalent of talking about u’v but computing the equivalent, v’u.

Why might we view it that way? Among other things, we can apply u to v without even having a dot product!

Let’s return to finding the dual basis E.
end digression

If A is square and of full rank — i.e. we really do have a basis, hence a linearly independent set of vectors — then we can invert A and write

.

In this case, however — with A our color matching functions –as we first saw in the first Cohen post, A is not square and therefore cannot possibly be invertible. The 3×3 matrix A’A (, however, is invertible, and we construct the pseudo-inverse

and substitute it for , getting

.

(We remember that A’A is symmetric, i.e. its own transpose.)

BTW, it is very convenient to write the equation that way, because to check that we have a dual basis we want to compute . We might as well have found first.

When we check it, we see that .

Let me call the three dual basis vectors E1, E2, E3, and let’s look at them. Here they are in columns…

… and here is what they look like (yellow, violet, orange in order):

If we multiply each of those spectra by our original numbers {356.817,354.641,271.109}… and add the results.. we get

Vdb = {0.240105,2.51498,23.8508,64.4948,61.7507,33.9437,25.2229,43.8813,
66.2108,82.129,90.1716,83.8477,60.51,30.3654,10.9919,3.10251,0.753495,
0.189327,0.0489062,0.0147546,0.}

which looks like this (you can see that the orange one shrank a little compared to the other two):

where I have shown the scaled basis vectors (i.e. the scalar component times each vector), and added their sum in white.

Let’s look quickly at the reflected spectrum (in yellow) and that sum (in white):

They are quite different. What’s going on?

What does the dual basis describe?

Well, look at the long wavelengths of the white curve. We — i.e. the standard observer — can’t see that part of the reflected spectrum at all. Parts of the reflected spectrum are invisible to us.

Duh. Been there, done that. “Metameric black”.

What we have — a basis of 3 vectors — is a basis for the visible-to-us part of the reflected spectrum — i.e. a basis for what we called the fundamental. We have in fact found that part of the reflected spectrum which the standard observer actually sees; and by subtraction we could find what part of it is– let’s call it the residual — what we have called a metameric black. (“Metameric” because it is a nonzero spectrum; it is perceptually black, but it is not physically the absence of light.)

We have the fundamental; it’s that white curve. But I haven’t really proved that.

It’s not hard. By definition, each of the 3 vectors in the dual basis is in the pre-image of . And, since has 3 rows, it is of rank at most 3, so its pre-image is of dimension at most 3. In fact, we know that this is of full rank, hence of rank 3, and therefore its pre-image is of dimension 3.

If the 3 dual basis vectors are linearly independent (i.e. the E matrix is of full rank 3), they span a 3D space which is a subspace of the 3D pre-image — so they span the entire pre-image.

And I’m certain that being of full rank guarantees that the dual basis is of full rank.

But exactly how do I prove that? Oh, the pre-image of is the range of A, hence of dimension 3.

Do you find it convincing?

Let’s look at it another way. But remember that we really did just get the fundamental f of the reflected spectrum S, and we could compute the residual n by subtraction; the residual is the reflected spectrum minus the fundamental.

n = S – f.

An orthonormal basis for fundamental spectra.

For the other way, we need a really nice basis for the entire domain of the matrix, one which decomposes its domain into nullspace and not-nullspace (i.e. nullspace and preimage of the range). But we know how to do that.

The SVD (Singular Value Decomposition).

We will get an orthonormal basis v1 for the pre-image of , then a projection operator R onto that pre-image, then get the fundamental by applying the projection operator to the reflected spectrum.

First, let

.

(The first time I did this, I let . Now we’re more experienced, but if you want to use B instead of , feel free.)

Now, v is 21×21. Then the leftmost 3 (= rank of ) columns of v are an orthonormal basis for the preimage of , and the rightmost 18 (= 21-3) are an orthonormal basis for the nullspace of .

Here are the leftmost 3 columns of v:

We should look at those. First alone:

… and then juxtaposed with the dual basis:

OK, the orthonormal basis is different from the dual basis. (Actually, one pair look similar. I’ll take a look in a moment.) These vectors (spectra) are orthonormal, the dual basis is not. (To check that, just compute ; it’s a 3×3 identity matrix.)

Let me look at the two basis vectors that were similar. From the dual basis, the first vector; from the orthonormal basis, the 3rd vector.

Interesting. I do not know that this result has any significance, but what’s the point of drawing pictures and not looking at them?

Let’s continue.

The projection operator R.

Having an orthonormal basis for the preimage of A^T, we know that we can construct a projection operator R from the domain of onto the preimage; we just take the product of v1 and in the other order:

It’s 21×21, so I’m not going to show it to you. But, for checking, I’ll show you the first column:

And I’ll tell you that the last column is identically 0.

The fundamental, again? and more.

So. Apply the projection operator R to the reflected spectrum S to get the fundamental f…

R S = f = {0.240105,2.51498,23.8508,64.4948,61.7507,33.9437,25.2229,43.8813,
66.2108,82.129,90.1716,83.8477,60.51,30.3654,10.9919,3.10251,0.753495,
0.189327,0.0489062,0.0147546,0.}

and subtract to get the residual n = S – f:

S – f = n = {4.85989,26.2994,19.1132,-14.6673,-7.32707,18.6749,27.0703,15.2259,
3.00644,-13.029,-23.7822,-12.7477,13.5965,41.0307,59.1247,63.6874,
60.2497,52.1091,59.2801,33.022,39.625}

We can check that n is in the nullspace of by computing :

We should also check that the fundamental has the same 3 values as the entire reflected spectrum. That is, we had

for applied to the reflected spectrum. If we now apply to the fundamental, we get

The same. Of course. (I hope you think so.)

Let’s look at the decomposition. The reflected spectrum is in yellow, the fundamental in white, and their difference n, the residual, is in black.

Yes, where the fundamental is zero, the reflected spectrum and the residual coincide.

At the risk of annoying you, let me emphasize: the standard observer (at 20 nm intervals!) would perceive the yellow spectrum and the white spectrum as being exactly the same color. We would say they are metameric spectra. He is actually perceiving the yellow spectrum, but if we somehow generated the white one, it would look the same to him.

And now we should compare the fundamental with the vector Vdb constructed from the dual basis. We could look at the numbers, of course, but an easy way to see that they coincide is to plot only points for Vdb and only a piecewise-linear curve for the fundamental.

They are the same.

Thus we confirm that the 3 values are the components of the fundamental (not of the entire reflected spectrum) wrt the basis E dual to A.

The dual basis tells us what those first 3 numbers mean, and would let us construct the fundamental, and thence the residual.

The orthonormal basis, on the other hand, let us compute a projection operator, which gets us the fundamental by applying the projection operator to the reflected spectrum.

If i could only do one of these, it would be to construct the dual basis rather than the projection operator… but computing power is cheap and there’s no reason to do only one of these.

Incidentally, we could have simply appled the Gram-Schmidt orthogonalization procedure to the dual basis… but I really like having a basis for the null space, and that’s what the rest of the v matrix from the SVD got me.

I’m going to wait, however, before looking at the basis vectors for the null space.

Color: from spectrum to tristimulus

Introduction

Edit 30 Oct 2009:
This is the first post about an example on p. 160 of Wyszecki & Stiles (see the bibliography).

and find another edit about the components wrt the dual basis.
End edit.

As always, if my chatter is confusing, just start following the examples.

In a sense, what I am about to do now is one of the most important color calculations I will show you.

How do we get from a spectrum to tristimulus values?

I am going to get from a spectrum to XYZ values, and then to xyz values. Recall that XYZ come from xyz bar tables, which in turn come from rbar, gbar, bbar tables by a change of basis. Then xyz values come from XYZ by making the sum = 1.

For more about the CIE, see the previous post.

The question before us is to get from a spectrum to XYZ (and then to xyz).

What we are talking about is the perception of color, rather than the reproduction. Sort of. We wil use “the standard observer” rather than my eyes or your eyes, but we are going to translate an arbitrary spectrum into 3 scalar values, which can then be transformed to some RGB values in a color space. But not in this post.

This example comes from w&s (Wyszecki & Stiles, in the bibliography), page 160. Believe it or not, I am going to work it four ways — but I will give you the full details for only two of them, pictures for all four of them.

The original problem uses data at 10 nm intervals, and it uses the 1964 CIE. Working this example as they did showed me that there were 3 typos in their solution.

First, however, I will work the problem using the data at 20 nm intervals, with the 1964 CIE. This will show you explicitly how the calculations proceed. At 20 nm intervals, I can display the data itself.

Second, I will work the original problem as they posed it. I will not get their answers, so I will also show you their typos.

Third and fourth, I will again work the problem at both 10 nm and 20 nm intervals, but using the 1931 CIE. As will be my ongoing custom, I will show you all the numbers only for the 20 nm case. Ultimately, what choice we make isn’t going to matter a lot. And that’s worth knowing.

We are given a reflectance spectrum for a plastic object illuminated by daylight. In fact, we will use “standard illuminant D65″ for “daylight”. The “reflectance spectrum” is the reflection from a perfectly uniform — perfectly equal energy — illuminant.

Let me say that another way. We are going to assume that the object is illuminated by one spectrum (“daylight” D65) and that the reflected spectrum which we perceive is the combination of the D65 spectrum and the so-called reflectance spectrum of the object.

The reflectance spectrum is a property of the object; by combining it with any arbitrary illumination spectrum, we get the reflected spectrum.

We are asked, in the given example, to find the XYZ tristimulus values, and the xy chromaticity coordinates, using the 1964 CIE color matching functions. That is, we will use the 1964 xyzbar tables.

(I didn’t plan to show this to you, but prudence dictated that I at least work it out. When it turned out that they had slightly wrong answers… well, I have to show it to you.)

Here’s what we will do. We will find the reflected spectrum, i.e. what is reflected from the object illuminated by D65 rather than by an equal energy illumninant. Then we will apply the xyzbar color matching functions to the reflected spectrum to get XYZ tristimulus values. Then we will get the normalized xyz values.

For more about the color matching functions, our matrices A for each problem, see the first Cohen post, and the second Cohen post.

CIE 1964 at 20 nm

Here is the reflectance spectrum , from 380 to 780 at 20 nm intervals. As will be my custom, I am taking a subset of their values — I am not averaging the numbers. That is, their example provides data every 10 nm, but I am using half of that given data. For one thing, as I said, that gives us matrices of manageable size.

Here’s a picture of that reflectance spectrum:

A word about scaling. It is appropriate that the reflectance spectrum be less than, or on the order of magnitude of, 1. At any given frequency, the amount reflected can range from significantly less to significantly more (since light absorbed at one frequency can be emitted at another) than the incident spectrum — but we dont want to be multiplying by a hundred. That would be bad.

Here is the D65 spectrum, from 380 to 780 nm at 20 nm intervals:

Another word about scaling. The D65 spectrum ranges from 50 to more than 100. I have no idea what units it is measured in — but the reflected spectrum will be in the same units, whatever they are.

Here is a picture of that illumninant spectrum:

As I said, the scale is 100 times that of the reflectance spectrum. This doesn’t matter for the computations. If, however, I want to show both of these on the same axes, I will need to scale these values — but for the graph only, not for the computations.

Here are the 1964 xyzbar tables from 380 to 780 nm, at 20 nm intervals. This is the A matrix for this problem.

And here is a picture of the color matching functions:

Those are our ingredients: the reflectance spectrum, the illiminant spectrum, and the color matching functions.

Now, let’s get the reflected spectrum. We multiply the values pointwise… i.e. at every wavelength, what is actually reflected is the product of the reflectance spectrum and the illuminant spectrum.

(Perhaps I should be calling this the perceived spectrum, since it is what our standard observer is seeing. I should not call it the incident spectrum, because that would be a sensible alternative to “the illimination spectrum”, in this case D65. Incident illumination on an object gives us a reflected spectrum, which our standard observer perceives.)

The result of multiplying the reflectance spectrum and the illuminant spectrum D65 is

and the 3 spectra look like this (after I divide the D65 and S by 100, otherwise the reflectance spectrum is effectively the x-axis):

Let’s be clear: the reflectance spectrum is in green, the illuminant spectrum is in blue, and the product, at each wavelength, is the reflected spectrum, in black.

Before we deal with the reflected spectrum, however, I want to deal with the illumination itself. If we apply the color functions to the D65 spectrum, that is we compute

,

we get

{554.645, 582.565, 631.147}.

Are those our tristimulus values? No way. These numbers depend on the scale of the spectrum, and the number of components of the vectors. Somehow we have to normalize them, and we will do that for the total energy received.

As an aside, which I will try to remember to repeat in subsequent posts…. We have spectra and color-matching functions which live in vector spaces. The result of applying color matching functions to spectra is a set of 3 scalars, i.e. 3 real numbers. But our vector spaces are not bounded spaces, so those 3 numbers are not, in principle, bounded. They are simply the components of a spectrum wrt a basis. Edit: oops. Careful! They are components of the fundamental of the spectrum. End edit. OTOH, we have color spaces, such as RGB, in which the RGB values are bounded (whether between 0 and 1, or 0 and 100, or 0 and 255, they are limited).

We need to move from the unbounded results of vector operators to the bounded values of color spaces.

But how do we do that?

Well. Recall that ybar, the second column of A, is the photopic response (the total response of the cones) of the standard observer. This is incredibly useful: the second column of AT tells us the total energy perceived. If we were doing calculus, we would compute the areas under all three curves. But we’re doing linear algebra, and the appropriate number is the middle one.

We get XYZ values for the illumnination perceived by the observer by dividing those 3 numbers by the middle number, 582.565:

XYZ = {0.952075, 1., 1.08339}.

Let’s keep going with this. How do we get xyz coordinates from XYZ?

Normalize the sum of the XYZ to 1. The sum is 3.03547 … so we get that x,y,z are

xyz = {0.31365, 0.329438, 0.356911}.

Since all three values are close to 1/3, this is a plausible number for illumination.

Now we can do the same thing with the reflected spectrum. That is, we apply the matrix to the reflected spectrum S:

{386.796, 381.417, 293.87}.

We are not going to divide by this middle number — this is the reflected spectrum, not the illimination spectrum. We need to divide by the same number as before, i.e. by the photopic response to the illumination, not the photopic response to the reflection. That number was 582.565, and now we divide our 3 values by that to get X,Y,Z…

XYZ = {0.663953, 0.654719, 0.504441}.

And, as before, we compute xyz values, too: we normalze the XYZ so that their sum is 1. Since the sum is 1.82311 … we get

xyz = {0.364186,0.359122,0.276692}.

This isn’t all that far from the illuminant.

Those dots are the outline of the 1964 chromaticity diagram. I got it the same way as the 1931, just with different data. Get the xyz tables — the ones where the rows add up to 1 — and just plot y versus x.

The red, green, and blue dots show the pure primaries for the 1964 standard. The white dot is (x,y) for D65, and the yellow dot is (x,y) for the reflected sprectrum. (Oh, I don’t think the reflection is yellow, but it’s rather closer to the light source than to much else.)

CIE 1964 at 10 nm

Now let me show you the pictures, and the answers, for 10 nm intervals. In fact, let me juxtapose the new pictures (on the right) with the old (on the left):

Here are the old and new reflectance spectra:

Here are the D65 spectra:

And here are the color matching functions (xyzbar tables) from 380 to 780 nm The data on the right is the A matrix for this problem.

We have the same ingredients: the reflectance spectrum, the illiminant spectrum, and the color matching functions.

Now, let’s get the reflected spectrum. We multiply the values pointwise… i.e. at every wavelength, what is actually reflected is the product of the reflectance spectrum and the illuminant spectrum.

The result of multiplying and D65 together looks like this (after I divide the D65 and S by 100, otherwise the reflectance spectrum is effectively the x-axis):

As before, the reflectance spectrum is in green, the illuminant spectrum is in blue, and the product, at each wavelength, is the reflected spectrum, in black.

First, we apply the color functions to the D65 spectrum, that is we compute

,

and we get

{1102.11, 1162.02, 1247.78}.

Here we see quite clearly that these numbers depend on the scale of the spectrum and the number of values: they are about double what the previous values were. Oh, of course: our vectors are twice as long. (Even though the XYZ and xyz values will be comparable, these initial 3 values will not always be comparable. As I said, they depend on the dimension of our spectra and on the scaling of the illuminant spectrum.)

Once again, we need to scale those numbers, dividing each of them by the middle number, which is the total perceived illumination.

So we get XYZ values for the illumninant by dividing those 3 numbers by the middle number, 1162.02:

XYZ = {0.948445, 1., 1.0738}.

How do we get xy coordinates? Normalize the sum to 1. The sum is 3.02225 … so we get

xyz = {0.313821, 0.330879, 0.3553}.

As before, this is a plausible number for illumination.

Now we do the same thing with the reflected spectrum. That is, we apply the matrix to the reflected spectrum S:

{769.263, 761.189, 581.056}.

Once again, we divide by the same number as before, i.e. by the photopic response to the illumination, not to the reflection. (That number was 1162.02.)

Now we divide our 3 values by that to get X,Y,Z…

XYZ = {0.662005, 0.655057, 0.50004}.

And, as before, we compute xy values, too: we normalze the XYZ so that their sum is 1. Since the sum is 1.8171 … we get

xyz = {0.364319, 0.360495, 0.275185}.

You can’t tell much from these two pictures; the final picture in this post will plot (almost) all of the xy answers. The numbers using 20 nm intervals were not much different: for the XYZ values of the reflected spectrum, we had:

XYZ = {0.663953, 0.654719, 0.504441}.

The corresponding answers using 10 nm intervals were:

XYZ = {0.662005, 0.655057, 0.50004}.

Pretty darn close. 10 nm or 20 nm, both using CIE 1964, not a whole lot of difference.

Now, what were their answers? They got sums of

{758.6, 760.8, 580.8}

… while I had

{769.263, 761.189, 581.056}.

We do not agree. What’s going on? I found that there is a problem in their work at 640 nm.

For a wavelength of 640 nm, we agree on the values of all the ingredients. We agree that the reflectance spectrum value is 0.853 … the D65 illuminant spectrum value is 83.7 … the values of that row of A are {0.4316, 0.1798, 0.}.

The pointwise products of these numbers are…

{30.8146, 12.837, 0.}

but what they show is

30.2 12.6 0.

Two typos. I have no idea why. And they appear to be the only typos pertaining to individual wavelengths. Of course, it is possible that the error is in the reflectance spectrum instead, but I have no other data-source to compare that to. In any event, it is clear that their entries (30.2 and 12.6) are not in fact the product of the inputs they show.

There’s another problem. I pasted the column that leads to the X-column-sum into a spreadsheet and computed the sum. For the data they have shown, the first sum is really 768.6 rather than 758.6 . That’s the third typo.

(In addition, they round off at each wavelength, before adding the columns: this leads to some round-off error. That’s why we differ on the sum that leads to Z; no typos there, they just rounded before adding.)

But, where they show X,Y,Z as (change of scale, decimals versus percentages; ignore it) …

XYZ = {65.3, 65.5, 50.}

the correct answers are

XYZ = {66.2005, 65.5057, 50.004}

and even using 20 nm intervals instead of 10, I got

XYZ = {66.3953, 65.4719, 50.4441}.

For x,y, they show

xy = {0.3612, 0.3623}

but correct (at 10 nm) is

xy = {0.364319, 0.360495}.

Just working at half the resolution (20 nm), I got

xy = {0.364186, 0.359122}.

I daresay their mistakes are worse than my cutting the resolution in half.

CIE 1931 10 nm

To switch from 1964 at 10 nm to 1931 at 10 nm requires only that I change the A matrix. (I have everything else in memory at 10 nm intervals. The reflectance, reflected, and D65 spectra do not change.)

Let me show 1931 on the right, for comparison with the 1964 color matching functions on the left. Note that the vertical scales are different.

Although they have not changed, I show again the reflectance spectrum and the illuminant spectrum for this case:

Go ahead, curse me out. After showing different cases side-by-side, I have just switched and shown side-by-side graphs for two different things for one and the same case.

Those are our ingredients: the reflectance spectrum, the illiminant spectrum, and the color matching functions. Only the color functions have changed, from 1964 to 1931.

Now, the pointwise product of the illuminant and reflectance spectra will be exactly the same reflected spectrum.

Let’s be clear: the reflectance spectrum is in green, the illumninant spectrum is in blue, and the product, at each wavelength, is the reflected spectrum — in black.

As before, we apply the color matching functions to the D65 spectrum, that is we compute

,

getting

{1004.46, 1056.89, 1150.02}.

These numbers are comparable to the previous case: we are using vectors of the same dimension.

Once again, we need to scale those numbers, dividing each of them by the middle number, which is the total perceived illumination.

So we get XYZ values or the illumninant by dividing those 3 numbers by the middle number:

XYZ = {0.950398, 1., 1.08812}.

To get xy coordinates, we normalize the sum to 1. The sum is 3.03852 … so we get

xyz = {0.312783, 0.329108, 0.358109}.

Now we do the same thing with the reflected spectrum. That is, we apply the matrix to the reflected spectrum S:

{707.724, 707.564, 537.101}.

Once again, we divide by the same number as before, i.e. by the photopic response to the illumination, not to the reflection. That number was 1056.89 .

So we divide our 3 values by that to get X,Y,Z…

XYZ = {0.669631, 0.669479, 0.508191}.

And, as before, we compute xy values, too: we normalze the XYZ so that their sum is 1. Since the sum is 1.8473 … we get xyz as

xyz = {0.362491, 0.362409, 0.275099}.

Again, this isn’t all that far from the illuminant. We’re getting very similar answers.

The CIE boundaries look pretty similar, but the triangles — the gamut of the 3 primaries — are different. And recall that the line x + y = 1 is tangent to the right-side boundary; furthermore, for the 1931 CIE, one side of the triangle is very, very close to the tangent line.

CIE 1931 20 nm

Having shown you pictures with the 1931 CIE and 10 nm intervals, let me close by showing you the numerical data, as well as pictures, for the 1931 CIE 20 nm solution.

To switch to 20 nm intervals, I need to rebuild the color matching functions. Here’s the data:

and here are the graphs (20 nm and 10 nm side by side):

I need to restore the reflectance spectrum and the illuminant spectrum. As you recall, they were

and

The reflected spectrum was the pointwise product of those, namely

Those are our ingredients: the reflectance spectrum, the illuminant spectrum, and the color matching functions. Only the color functions have changed, from 1964 to 1931.

The product of illumninant spectrum and reflectance spectrum is the same as shown before for 20 nm intervals:

As before, the reflectance spectrum is in green, the illumninant spectrum is in blue, and the product, at each wavelength, is the reflected spectrum — in black.

As before, we apply the color matching functions to the D65 spectrum, that is we compute

,

getting

{506.693, 529.764, 581.089}.

These numbers are smaller than the previous case: we are using vectors of half the dimension.

Once again, we need to scale those numbers, dividing each of them by their middle number, which is the total perceived illumination.

So we get XYZ values or the illumninant by dividing those 3 numbers by the middle number, getting:

XYZ = {0.956451, 1., 1.09688}.

We get xy coordinates by normalizing the sum to 1. The sum is 3.05333 … so we get

xyz = {0.313248, 0.327511, 0.359241}.

Now we do the same thing with the reflected spectrum. That is, we apply the matrix to the reflected spectrum S, getting:

{356.817,354.641,271.109}.

Once again, we divide by the same number as before, i.e. by the photopic response to the illumination, not to the reflection. That number was 529.764, so we get XYZ are

XYZ = {0.673539, 0.669432, 0.511754}.

And, as before, we normalze the XYZ so that their sum is 1. Since the sum is 1.85472 … we get

xyz = {0.363148 ,0.360933, 0.275919}.

OK, let’s compare.

Here are the unscaled sums (including their incorrect answers). Reading down, these are my 1964 at 20, then at 10, their answers, my 1931 at 10, then at 20.

and the XYZ values scaled by the photopic response…

and the xyz values

There’s not a whole lot of difference there, whether we use 1931 or 1964, or 10 nm or 20 nm, or make a few mistakes!

To look at it, here are my 4 computed xy, “white points”, for the illuminant, and the 5 xy points for the reflections.

(I didn’t grab their white point — the x,y coordinates — for D65, so only my 4 computed values are shown.)

The most important single computation in all this was to scale the values by the middle sum, i.e. by the total energy perceived by the observer. To be specific, we applied the ybar function to the illuminant spectrum.

I know, it looks just like the dot product of two vectors. But Cohen taught us that the color functions are really linear functionals. We want to visualize that computation as the matrix product of a row (the ybar) and a column (the D65 spectrum).

Of course, only such a trick as that division by the middle number could distract us from the marvelous computation of 3 numbers in the first place. We get from a spectrum to 3 scalars by applying 3 color matching functions (the transpose of our A matrix) to the reflected spectrum, which is also the spectrum incident on the eyes of the standard observer.

The first 3 numbers we get are, as we saw before and will see again, the (potentially unbounded) components of the incident spectrum wrt the dual basis spectra. The second set of 3 numbers are tristimulus values in XYZ space. The third set of numbers were xyz coordinates for the chromaticity chart.

The first 3 numbers depend on the dimensionality of our spectra, how many samples we have. In fact, we compute two such sets: one for the illuminant and one for the reflection. The middle number of the 3 for the illuminant sets the scale for white = (1,1,1).

The second set, XYZ, is found by dividing the first 3 numbers by that middle number. That is true whether we have XYZ for the illuminant or for the reflection: divide either by the total energy of the illuminant.

The third set, xyz, is found by dividing the XYZ numbers by their sum — instead of by their middle, or anybody’s middle.

I must emphasize that we have seen 4 different choices for the A matrix, and those were all “xyz bar” values! “The” A matrix is hardly unique.

More to come. You better believe it.

Color: the 1931 CIE color-matching functions and chromaticity chart

Introduction

The CIE chromaticy chart is one of the things we are headed for. Here is a black-and-white drawing of its boundary. I’ll show you later how I got it.

You should be able to find color versions of the CIE chromaticity chart (1931) all over the Internet. Here’s one that shows some curves of constant hue. I’ll point out now that these curves are not straight lines. I will also use that image to decide that something is a greenish-yellow or yellow-green.

Wikipedia seems to have a lot of information, such as this, for example.

What my drawing shows is that I can reproduce the analytical content of it. Black and white, but no color.

How did I do that?

the color-matching functions and chromaticity coordinates

I have four tables of data, from Wyszecki and Stiles, for which see the bibliography; and I will abbreviate them “w&s” in this post.

“rgb” are chromaticity coordinates (row sums = 1), “rgb bar” are color-matching functions. The full heading on the combined “rgb” and “rgb bar” tables is “Chromaticity Coordinates and Color-matching Functions of the CIE 1931 Standard Colorimetric System, with respect to real primary stimuli:
R at nm;
G at nm;
B at nm”

We also have “xyz” chromaticity coordinates (row sums = 1), and “xyz bar” color-matching functions.

All 4 tables are at 5 nm intervals from 380 to 780.

Well, these 4 tables come in a variety of sets. I chose to use the 1931 tables at 5 nm intervals.

They have been interpolated to 1 nm intervals. They have been corrected by several authors in different ways. They were re-done in 1964 by the CIE to use a 10° field instead of the 2° field of the 1931 standard.

My choice of the uncorrected 5 nm 1931 standard is not universal — but it is an extremely common choice, and an awful lot of the work done since 1931 has used it.

At my level of usage, it probably doesn’t matter a whole lot. In fact, a whole lot of other things will turn out to be more important.

Since I’ve had more than a little time to get used to the 4 different tables, and you may not have, let me list them:

• rgb bar
• rgb
• xyz bar
• xyz

My source for the tables, as I said, was w&s. I am not going to reproduce them — as far as I can see the CIE charges for some of them — but I wouldn’t be surprised if you can find them on the Internet.

In principle, three of the tables are secondary. Their source is the table of rgb bar. First, let me show you what they look like.

Here is a picture of the rbar, gbar, bbar data. The x-axis runs through wavelengths from 380-780, the points are 5 nm apart.

It is clear that rbar is significantly negative over part of its range (440-545 nm); it is true, but not clear, that gbar is very slightly negative over 380-435, and bbar is very slightly negative over 550-655 — and it’s shown as 0 for longer wavelengths.

In fact, at any given wavelength, exactly one of the three is negative; nevertheless, only the red rbar gets seriously negative.

What do these negative numbers mean?

I suspect that the simple description of the color matching experiment is too simple experimentally, but I also believe it suffices in principle. Imagine that you are looking at a split screen, at two halves of a colored disk. The left half of the disk is a monochromatic color. (Think of a tunable laser — but of course they didn’t have lasers in the 1920s when the experiments were done!) The right half of the disk is controlled by the viewer, who has three dials and may choose mixtures of red, green, and blue in order to make the two halves of the disk match.

It turns out that we can’t do that. Over most of the range, we have to change the color of the left half of the disk as well as the right half. We can’t match a wavelength of 380 nm (the violet end of the visble spectrum) with red, green, and blue. But if we add a tiny of green to the violet, we can match that combination with red and blue.

We record a negative number for green at that wavelength.

Until the blue goes to zero thru negative values, exactly one of the 3 numbers is negative. At long wavelengths, when blue is zero, we just use red and green.

So much for the negative values.

Oh, having said that I will not reproduce the table, I will present — so you can do these calculations yourself — one quarter of it. Here is the rgb bar table at 20-nm intervals.

The second table, rgb, is derived from the first in principle — but not in practice!

In principle, take the rgb bar table and recsale it by setting every row sum to 1.

In parctice, we don’t have enough digits to do that. They provide a table. It is called the rgb (chromaticity coordinates).

Here is what the rgb table looks like at 5 nm intervals.

And here is a subset of the rgb table, at 20-nm intervals.

Let me show you an example of why we need their table, instead of trying to generate it ourselves. At 380 nm, the rgb bar table has values…

which sum to 0.00119 .

If I divide that row by its sum, and round, I get…

but the corresponding row of the rgb table is…

.

Close, but not the same.

That’s two of the four tables. As it happens, they’re the two we may be able to dispense with in practice — but in principle, rgb bar is the starting point.

Let me emphasize that these values depend on a specific choice of red, green, and blue.

We can fix that.

The CIE defined a specific change-of-basis. The tri-stimulus coordinates (rbar, gbar, bbar) or (R,G,B) are a vector in R^3; they defined a vector denoted either (xbar, ybar, zbar) or (X,Y,Z).

Oh, although the table is denoted rgb bar, specific coordinates wrt it are usually denoted RGB; similarly, specific coordinates wrt the xyz bar table are usually denoted upper-case XYZ.

Bear in mind that these RGB are not the same as coordinates in RGB-color space. (The latter is bounded by 0 and 1; the former is not!) We’ll see more about this down the road.

Let me show it to you:

And let me show you the cut-down data, every 20-nm; this is a subset of the “xyz bar” table:

Here is the change-of-basis matrix:

They tell us we can check it, for example, by applying their matrix T to RGB coordinates, such as for example at 500 nm…

Apply T, and round off, getting

… and compare that to the XYZ coordinates at 500 nm:

.

The same.

Still… raise your hands and step away from the keyboard.

What we just did was apply their matrix T to RGB components to get XYZ components. Perfectly natural. But we know that the transition matrix delivers old components when it is applied to new ones. If we want to call that given matrix T the transition matrix, then we must think of XYZ as the old components.

In practice, this really just means that we better know which way the transformation goes. They gave us the transition matrix for old components XYZ and new components RGB — even though they derived the XYZ basis from the RGB basis.

OK? If not, please be patient.

Now just where did that change-of-basis transformation come from? Malacara has a readable discussion. The most important point, to my mind, is that the Y column vector, the second column, is chosen to be the “phototopic response” (1924) of the human eye, that is, daylight vision, that is, color vision. There’s not much point in my digging out that data and plotting it: it is exactly the green curve in the drawing of the xyz bar table.

Second, all the values X,Y,Z are non-negative. Third, they also set the column sums almost equal to the value for the Y column. There was another criterion, but I’ll show it to you after I show you the fourth table.

Table 4 comes from table 3 as table 2 comes from table 1: get table 4 by setting the row sums of table 3 to 1, just as we got table 2 by setting row sums of table 1 to 1. All in principle, of course, or having more digits available than were printed.

Here is what the xyz table looks like:

and here is the cut-down data:

Just for the fun of it, check the column sums of XYZ. I said they were almost equal. I get

.

I do not know if that is deliberate, or whether its round-off and the real tables have more digits. (I don’t really care. But I couldn’t very well say they were equal….)

The CIE chromaticity chart

Ok, here’s the next major idea. Since x+y+z = 1, this table is essentially 2D, and we can plot y vs x.

That is exactly the picture I started this post with, the boundary of the chromaticity chart.

That’s all it is: take the first two columns of the xyz table, read each row as (x,y) coordinates, and plot each row.

The fourth criterion was to make the line x+y = 1 tangent to the curve. (In fact, the earlier criterion about X,Y,Z being non-negative was implemented by setting the y-axis tangent to the boundary.

Let me show you:

So, 4 tables, one 2D chart. We could also draw a 2D chart of g vs r, or b vs r, from the rgb table, since its rows each sum to 1.

What we have lost, in that chart, is the intensity. To recover XYZ from xy, we would need, for example, xyY; x,y and any one of X, Y, or Z would suffice.

Let me start using that chart right now. We have RGB basis vectors (1,0,0), (0,1,0), and (0,0,1}. What are their (x,y) coordinates?

There are a few ways to do this. Given the scale at which we are working, the easiest is to interpolate the x,y coordinates from the xyz table at wavelengths 435.8, 546.1, and 700 (that last one is a table value).

I get x,y coordinates of

R = {0.7347,0.2653}

G = {0.2737,0.717411}

B = {0.166541,0.00892768}

Let’ s plot them on the chromaticity chart… and let’s draw lines connecting them.

Think back to simplices. All of the points inside (and on) that triangle can be written as linear combinations of the coordinates of the 3 vertices; most importantly, the coefficients of those linear combinations are between 0 and 1, inclusive. That is, when the coefficients are in the closed interval [0,1], the combination is inside or on the triangle.

That restriction to [0,1] is what restricts us to the inside (and on) the triangle. That triangle is said to describe the gamut of these three R,G,B vertices.

The points outside the triangle have either negative coefficients (or perhaps coefficients larger than 1; I’m not sure). This is why we have negative values inside the primary table, rgb bar: all the points on the curved boundary are outside, or just on, the rgb bar triangle.

That triangle is said to to show the gamut of those choices for R,G,B primaries. One of the most common uses of the chromaticity chart is to compare gamuts.

From Foley, van Dam et al., I get the following x,y coordinates for the “standard NTSC RGB phosphor” (the color TV standard):

R = {0.67,0.33},

G = {0.21,0.71},

B = {0.14,0.08}

Let’s add them to the CIE, triangle in red:

A TV with these phosphors can reproduce colors in the red triangle.

Another set of coordinates from Foley & van Dam et al., for “short-persistence phosphors” are:

R = {0.61,0.35},

G = {0.29,0.59},

B = {0.15,0.063}

Let’s add them to the drawing, with their triangle in green:

Whatever these are, they can only produce the comparatively limited set of colors bounded by the green triangle.

To reproduce the entire CIE chart, within a triangle — i.e. with linear combinations having components in [0,1], we would need to have B at x,y = (0,0), R at x,y = (1,0), and G at x,y = (1,0).

Easy enough to draw, but don’t ask me how to get such sources physically.

A different form of the same principle — linear combinations of vectors — leads to the idea of mixing two colors. That is, instead of looking at all that can be produced by a particular set of RGB phosphors, we take two points and consider the line joining them.

A particularly useful choice is when one of the points is defined as “white”. Here are the x,y coordinates of a white known as D65, a standard form of “daylight”; here also are x,y coordinates of a form of green (from Andrew Glassner’s “Principles of Digital Image Synthesis”, volume 2, Appendix G; Morgan kaufmann, 1995), known as “foliage” on the MacBeth Color Checker.

p65 = {0.31271,0.32902}

pf = {0.4002,0.3504}

I want the equation of the line thru those points, and considering triangles, I’m going to write the equation as

line[t_] := t pf + (1 – t) p65.

Like our discussion of triangles, if t is restricted to [0,1], then we get the line segment bounded by the two points. Only if we let t range outside [0,1] do we get line segments that extend beyond the two points.

For what it’s worth, the numerical equation of the line is

.

Here are the two points and the line segment joining them:

But I’d rather do something else: extend the line to the boundary, and consider the green point to be a combination of the white point and the boundary point.

According to my poster of the 1931 CIE, a wavelength of 570 is on the border between greenish-yellow and yellowish-green. Go to the first link in this post. So our “foliage” at the red point is a mixture of white and a wavelength very close to 570 nm. (Or, Malacara has a version of this drawing.)

(I meant to show the two points as black and red in both drawings, but I forgot. I didn’t want to imply that the dot was the “right” color.)

One final point. We could go from HSB to RGB (for some standard set of RGB) and compute curves of constant hue. They are not straight lines. Rather than prove this by computation, let me just point you at the very first link again: the boundaries which it shows are not straight lines, and those boundaries certainly ought to be curves of constant hue, namely the boundary hue. (It ought to be easy enough to show this, but not now.)

It is crucial that the XYZ values are independent of any particular RGB values. Of course, they originally came from the spectrallly pure ones — but once we have the XYZ basis, we can go to any other RGB basis we choose.

Having seen other gamuts, we should expect that people get XYZ values from a spectrum. I’ll show you that next.

Color: re-doing Cohen’s example

Cohen’s example again

Let me now show you how I would do Cohen’s example. (His computations, pretty much, were the previous post.) I cannot over-emphasize that he deserves a lot of credit for getting the mathematics right, even if he didn’t name it correctly or do it beautifully.

I start with the A matrix:

This, as I have said, is a toy example of real color matching functions. Each column corresponds to one of the tri-stimulus values; each row represents a wavelength.

My summary of Cohen began with a different matrix — the E matrix — only because it was easier to enter into the computer. A is primary.

The next issue we have to cope with is that we really care about the transpose rather than about A. A’ will be applied to a spectrum, and will deliver three tri-stimulus values for it. I’m not convinced that any notational choice I make now will be without problems.

(As I said in the previous post, I use A’ for the transpose whenever it is convenient.)

I am going to let

.

B represents the linear operator I wish to study. Thus, B is

It maps to . At most, it can be of rank 3; and then it would have a four dimensional null space.

What shall we do with it? All together now, class: find its singular value decomposition (SVD).

,

where u and v are orthogonal (hence square) matrices, and w is as close to diagonal as it can be, given that it is the same shape as B.

The Mathematica command is

{u,w,v}=SingularValueDecomposition[B];

where the final semi-colon suppresses printing.

It is convenient to look first at the singular values themselves, the nonzero “diagonal” elements of the matrix w:

.

That all three are “significantly” larger than zero tells us that the matrix B (and therefore its transpose, A) is, indeed, of rank 3.

In case this is new to you, the w matrix is

(And if this is new to you, check the tag cloud and the categories. There’s a lot examples of the Singular Value Decomposition. I could almost rename this blog “Applications of the SVD”.)

What about v?

It is 7×7. We know from the discussion of the four fundamental subspaces (look at the blue text) that the 4 rightmost columns of v are an orthonormal basis for the null space of B, and the 3 leftmost columns are an orthonormal basis for the pre-image of the range of B.

That is, writing a partitioned matrix [v] = [v1 v2], we have

We should confirm that B.v2 = 0, i.e. every column of v2 is in the null space of B:

Once having a basis, v1, for the pre-image of the range, we know that we can construct a projection operator as v1 v1′:

Because I didn’t call it R, I can recall Cohen’s R matrix, computed using a pseudo-inverse. I assure you that it is the same, we have P = R.

On the other side, since B is of full rank, its range is equal to its codomain, both being , so the entire matrix u is an orthonormal basis for the range of B.

Now the thing to do is, as before, split some given spectrum into its fundamental and its residual (metameric black). We only really have one good test, an equal energy spectrum.

We get the fundamental of our test spectrum…

…by applying our projection operator to it:

and that is exactly what we got before (come on, P = R; we have the same projection operator.)

Then we get the residual n as the difference…

and that, too, is the same as before.

As a refesher, let’s compute the components of the fundamental wrt the v1 basis. Since v1 is orthonormal, we may find the components of the fundamental by taking dot products, the fundamental dotted into each basis vector in v1. But that’s just

which says that f can be written as the sum

of the following three multiples of the basis vectors– I’m going to round off for clarity –

.

The sum does indeed (almost) match our fundamental. (The discrepancy comes from my rounding off.)

.

The last thing we need to get is the dual basis E for the rows of B. I just go look it up in the previous post…

:

Check that B E = I; it does.

And why do we care about the dual basis, i.e. the columns of E? Well, if we apply B to our test spectrum, we get

.

What are those three numbers? They are components of our test spectrum wrt the dual basis E!

That is, the sum of the following three vectors (which are multiples of the basis vectors in E)…

is

… which, once again, differs from our fundamental

only because of the rounding used to make the numbers readable.

To say that again: the dual basis E gives us a geometric interpretation of the tri-stimulus values we computed by applying B (i.e. A’) to a spectrum. The tri-stimulus values are components of the spectrum (equivalently, of its fundamental) wrt the E basis.

So. What did Cohen and I have that was the same? The matrices A and E, and the projection operator R. But note that where he computed R as

R = E’A,

and that’s short for ,

I would use

.

I’ll admit that’s not much of a savings if we have both A and E.)

What did we have that was related? Where he had two orthonormal bases F1 and F2, I have one, v1. All three of these bases span the same subspace of .

Digression: ways to show that F1, F2, and v1 span the same subspace of R7

One way to show that is to show that the columns of F1 and F2 are preserved by the projection operator cobstructed from v1. We have

and we compute

And yes, their difference is the zero matrix.

Similarly for F2.

Alternatively, F1 and F2 should generate the same projection operator as v1 did.

I assure you that is the same as P and R. And similarly for F2.

Alternatively, again, we could simply write the columns of F1 (then F2) as a linear combination of the columns of v1, and solve for the components.

Alternatively, yet again, we could try to find a transition matrix relating F1 and v1, then F2 and v2. That’s awkward because they’re not square matrices. Guess what? Use a pseudo-inverse! (Very likely the next math post.)

I have no idea what is special about F1 and F2. They are orthonormal bases, like v1, for the pre-image of the range of B. Since they span the same space, they are related by rotations (and possibly a reflection).

end digression, return to main thread

Finally, what did I do that Cohen did not? The singular value decomposition, which got me a basis, v2, for the null space of B — which he did not have. And I daresay that v1 was easier to get than either F1 or F2.

I think having v2 is incredibly useful. I get a basis for metameric blacks.

(Of course, Cohen had computed the matrices Ma and Me, and Ga and Ge, but I don’t need them because I wasn’t computing F1 and F2.)

Let me summarize that.

Cohen

• A and E;
• Ma = A’A, Me = E’E, where Me = Ma^-1;
• Ga and Ge, Cholesky decompositions of Ma and Me;
• orthonormal bases F1 = A Ge and F2 = E Ga;
• the projection operator .

Rip

• A;
• A’ = B = u w v’;
• [v] = [v1] [v2] where both v1 and v2 are orthonormal bases;
• the projection operator R = v1 v1′;
• the dual basis E.

(Strictly speaking, A is the dual basis to E, but the relationship is symmetric. Nevertheless, the rows of A’ live in the dual space V*, while the columns of E live in the original space V.)

I dispensed with everything leading to F1 and F2. I use v1 in place of F1 and/or F2. I get a basis, v2, for the null space.

I think the key point is using the projection operator to split a spectrum into fundamental and null vectors. It’s a unique orthogonal decomposition into a vector which has the same tri-stimulus values as the original spectrum, and a vector in the null space (metameric black, with tri-stimulus values = {0,0,0}.)

(The decomposition is unique to A, but A itself is not unique. More to follow.)

Cohen: “Visual Color and Color Mixture”

Introduction

Oct 10 edit: the third heading has been changed to “Computing E from A”

I want to work an example from Cohen, “Visual Color and Color Mixture” (see the bibliography, and this “books added” post). I am not, however, going to do it exactly the way he did. Nevertheless, I will show you everything he calculated.

Because I want to get everything of his into one post, I will break this into small sections. I expect that my next post will show you what I would have done instead.

All of his matrices can be found on p. 70 of his book.

What we have here is an extremely small example to illustrate “color matching functions” applied to light spectra, resulting in three real numbers which we call R,G,B or X,Y,Z. This is a prelude to using real color matching functions on real spectra. I will refer to the 3 numbers I get during the course of this example as “R,G,B”, in quotes because this is a toy example, and because down the road I’ll be computing “XYZ tristimulus values” in preference to RGB.

The A and E matrices: computing A from E

He has two matrices denoted A and E. They are closely related, and the typing is simpler if I start with the E matrix:

Let me emphasize that the A matrix is primary, and we should have computed the E matrix from it.

Let me answer a question at the outset. What is the point of the A matrix? Applied to a spectrum — careful! it’s actually the transpose A’ that is applied to a spectrum — it generates R,G,B or X,Y,Z values. This will get us from light spectra to colors! And the matrix E contains spectra; more to follow about that, I assure you.

Given E, he would compute the A matrix as

where I will frequently use A’ for the transpose , simply because it is sometimes typographically more convenient — and it is an extremely common notation. But note that I will write for the inverse transpose. In general, I will use whichever notation is more convenient for any particular expression, and I may even mix them shamelessly.

We get (with 4 places, as he showed it)

Hey what? What did he accomplish with that? First, we may say that he found an A such that both

A’E = I,

and E’A = I,

i.e.

A straightforward way of stating those relationships is: A’ is a left inverse of E; A is a right inverse of E’. We should also observe that by transposing

A’ E = I

we get

E’ A = I,

i.e. the second of those relationships. They are equivalent.

Now, we computed A from E as

,

Please note that the existence of depends on E having more rows than columns, and on E being of full rank. The matrix E is of rank 3, and the 3×3 matrix E’E is invertible. The 7×7 matrix EE’ is not, and can never be (not with E having 3 columns). And, of course, the matrices A and E themselves cannot possibly be invertible, since they are not square.

The A and E matrices: computing E from A

Had we been given the matrix A, we might guess that we could have computed E from A by some similar formula. Let’s try the obvious,

and see what we get. We have

.

That doesn’t look very good — but wait a minute. We have spent a lot of time computing “dispersion matrices”, things of the form E’E — and they are symmetric. That is,

.

So we continue, getting

.

That’s what we were hoping for, that E is computed from A exactly the same way as A was computed from E. Given the matrix A we would compute the matrix E as

.

By the way, we saw, in that process, that

.

This matters later when we define Ma and Me. That latest equation will tell us that Ma and Me are inverses.

The “real” relationship between A and E: dual (reciprocal) bases

All of Cohen’s work in this example treats A and E as matrices. There is another interpretation, a crucial interpretation. I will return to this in another post, but I have to tell you what he “really” did.

We have seen this computation before. Some of you may be sick of it. Well, it isn’t quite what we have seen before.

What we have seen before is a square matrix — specifically a transition matrix defining a basis — and we have wanted to compute the reciprocal basis. I even emphasized that I do not remember the answer, but I remember how to get it. If P is the given transition matrix, and R is the transition matrix for the reciprocal basis, then I want to compute all the dot products of the columns of R with the columns of P, and I want the answers to be an identity matrix. That is, I want an identity matrix to be the matrix product of rows of R’ with columns of P,

R’P = I

so R is the inverse transpose of P:

.

Wait just a minute. For a reciprocal basis, we just wrote

R’P = I.

For A and E, we had

E’A = I.

Guess what? Each column of E is a reciprocal basis vector for the columns of A. What we have are three columns of A, 7-dimensional but spanning a 3-D subspace. And the three columns of E are a reciprocal basis.

For reasons that will be far more obvious later, here is a place where I would choose to call the columns of E a dual basis rather than a reciprocal basis. Computationally — in finite dimensional spaces — there is no difference, but conceptually it will help a lot even there. In fact, I would say that we need to distinquish them conceptually precisely because we can’t tell the difference numerically.

It also turns out that it’s better to think of the columns of E (column vectors) as the original basis, and the rows of A’ (row vectors) as the dual basis. It might be best to think of the pair, columns of E and row of A’, as a set of dual bases, plural. Each row of A’ represents a so-called linear functional, a linear operator mapping a spectrum (a column vector) to a real number (the suffix -al in functional says the output is 1D). If this is new to you, relax. We’ll see it again, with real matrices A and E and real spectra.

I have barely discussed the dual space (that’s where dual vectors live, while reciprocal vectors live in the same space as the original basis) in these posts. We will get to it. If you know about it, good: for color theory, think dual basis rather than reciprocal basis. If you don’t know about the dual space, it’s okay to think reciprocal basis, for a while.

Okay, I’ll tell you that if you’ve ever seen differential forms, you’ve seen dual bases:

I should remark that Cohen never uses either term, reciprocal basis or dual basis. He knows that the relationship between A and E is crucial, he just doesn’t use the mathematical name for the relationship.

Now is a good time to remark, as I have before for other authors, that Cohen was almost certainly breaking new ground. I will not fault him for not saying it “correctly”; he deserves far too much credit for getting the mathematics right in the first place.

So, we still have exactly what we started with, the A and E matrices. In fact, we start with either one of them and get the other.

His dispersion matrices Me and Ma, and their Cholesky (LU) decompositions

I am going to show you the rest of what he did. We will later do almost all of this differently, but I had to know what he was doing. And if you are reading Cohen, then I have to show you a better way to do what he did. He computed (and named) the two intermediate dispersion matrices:

Me = E’E

Ma = A’A

and then he focuses on the inverses and . Oh, he also knows that Me and Ma are themselves inverses, as I noted earlier when we worked out E in terms of A:

.

(That is, there are only two distinct matrices among these 4 names.)

Anyway, here are Me, and its inverse…

and, for example, Ma…

which is, as it should be, .

Although he did not describe it this way, he next computed two matrices by doing Cholesky decompositions of Ma and Me. I think it is safe to describe the Cholseky decomposition as a special case of the LU decomposition. It is the result of an LU decomposition applied to a positive define symmetric (Hermitian in general) matrix, and the resulting L and U matrices are transposes, so there’s really only one of L and U to be found.

If you are following along with Cohen, let me warn you that I am going to change notation. Why? Because there’s a little magic going on here and a slightly different notation will focus on it. Besides, his notation seems silly. He has A and E, ends up with F1 and F2 — and along the way he uses G and Gbar.

I’ll keep his F1 and F2, since they are among his final results, but I am going to use Ga and Ge for the intermediate G’s. Trust me: you want me to do that.

Let me just do them. For Ma, we get…

What he has found is a lower triangular “square root” of Ma. That is, we have Ga such that

Ga Ga’ = Ma.

(This is his Gbar.)

We do it for Me, too…

As before, we have gotten Ge such that

Ge Ge’ = Me.

Note that I have used Ga and Ge to denote what we took the Cholesky decomposition of: Ga came from A’A, and Ge came from E’E, which is the inverse of A’A.

His orthonormal bases F1 and F2

He then defines

F1 = A Ge

F2 = E Ga.

That is what I want you to see: he uses Ge with A, and Ga with E. This guarantees that F1 and F2 are orthonormal matrices. (See below.) And if I view Ge and Ga as transition matrices, then F1 is apparently a basis for the column space of A, F2 for the column space of E. I have to confess that I don’t really care about F1 and F2 — I can do far better.

But, here are Cohen’s F1 and F2:

Okay, just why did F1 and F2 turn out to be orthonormal? That is, they each satisfy F’ F = I:

Let’s see why. Let’s compute F1′ F1. We have

F1 = A Ge.

Then

F1′ F1 = Ge’ A’A Ge = Ge’ Ma Ge

.

It works out because Ma and Me are inverses, and Ge is used with A to get F1. A similar calculation works for F2.

I want to point out that F1 and F2 are not dual to each other; we have, for example, F1′ F2 != I:

Of course not. Just as an orthonormal basis is its own reciprocal basis, so it “is” its own dual basis, at least in a finite dimensional space. I waffle and use quoted “is”, because the two bases which are dual to each other live in different vector spaces — but if one is orthonormal, the dual has the same numerical components. We will see this explicitly. In this case, it means that F1 “is” its own dual basis, and F2 “is” its own dual basis.

Let me say that again, approximately but simply: F1 is dual to itself, and F2 is dual to itself; they are not dual to each other.

But today we know how to get such bases without the Cholesky decomposition. We also know that we can use the same orthonormal basis for both of those spaces.

I will do that in the next post.

Using the A matrix: applied to the columns of E

I ask again, what is the point of the A matrix? Applied to a spectrum — careful! it’s actually the transpose A’ that is applied to a spectrum — it generates “R,G,B” values. Of course, this A matrix is a toy, but better to start with 3×7 than 3×81 or even larger. We’ll get there.

There are two interesting possible spectra to which we might apply A’, even from a toy matrix A. The obvious possibility, to me at least, it to apply A’ to E. (That is, we are applying A to 3 spectra at once.) Now, we already know the answer, because E was chosen as the dual basis. That is, we have already seen that

A’E = I.

That is, if we let the columns of E be E1, E2, E3, then we have

A’ E1 = (1,0,0)

A’ E2 = (0,1,0)

A’ E3 = (0,0,1).

In words, E1 is a spectrum which leads pure red; E2 generates pure green; and E3 generates pure blue. More precisely, E1 (etc.) is a spectrum which generates our chosen red. I’ll make more sense of that when we have real A matrices. (Oh, yes, there are more than one.)

Using the A matrix: applied to an equal energy (constant) spectrum

The second possibility is an equal energy spectrum:

Let’s not worry too much about that negative value; what it means is that we cannot actually match the equal energy spectrum (with whatever our toy “R,G,B” light sources are); instead, we have to add red light to the equal energy spectrum, and match that resulting mixture.

Perhaps now is a good time to say that the key point here is that we get 3 values; our perception of color vision is specified by 3 numbers. In practice, “CIE tristimulus values X,Y,Z” will be our objective.

The projection matrix R

There is one final crucial question. A’ maps from 7D to 3D. It is of full rank 3, so it has a nullspace of dimension 4. (If you need, now would be a good time time to look at fundamental subspaces.)

There are a whole lot of spectra out there whose R,G,B values are {0,0,0}. More importantly, any given spectrum can be split into two vectors, one in the nullspace of A’ and the other in the preimage of its range (i.e. in the range of A).

How would we find that decomposition?

Any direct sum decomposition has projection operators associated with it, so let’s find the projection operator onto the preimage of the range.

Malinowski did that. I did it differently. Cohen does it the same way Malinowski did, but it looks more compact because he defines both E and A. Cohen writes

R = E A’.

If we use

to expand Cohen’s equation, we get

.

We have also seen that this is the “hat matrix”, which projects an observation onto the subspace spanned by the least-squares fit. That is, we have

,

and

,

so

and if we define the hat matrix H by

,

we have

yhat = H y.

The key is that H is to X as R is to A.

Just for completeness, I compute R as Cohen did:

Since we know that projection operators are idempotent, we can check that . (It is.)

The projection of the columns of E

Having gotten the projection operator R, let’s use it. First, apply it to E, the dual basis to the rows of A’.

I hope it is not a surprise that we just got E back: the projection of E is E.

But what does that mean? No part of any column of E is in the nullspace of A’.

The projection of the equal energy spectrum

Recall the equal energy spectrum:

I will denote the projection as f, for fundamental.

Well, that has a nontrivial component in the nullspace, call it n. (We know that, as soon as we see f != ee.) We have

ee = f + n,

where n is simply the difference, n = ee – f, namely

Direct computation confirms that n is in the nullspace of A’:

Hang on just one second. What would you call a color whose “R,G,B” values were all zero?

I would call it black.

Because n itself, however, is a nontrivial spectrum — in the sense that there’s light there at each of those 7 frequencies! — let us refer to that spectrum as metameric black. (Two different spectra are said to be metamers or metameric if they appear to be the same color under the same illumination.)

It’s black to our eyes. (Although I could guess that if it were sufficiently intense, we might find ourselves with sore eyes from the radiation that is hitting them.)

I would emphasize that n depends more on A’, in a sense, than on the original equal energy signal ee. Of course, it came from ee and its fundamental, but once I have this n, I could add it to any other spectrum (of length 7!) without affecting the observed color of that spectrum. (In fact, it will depend on the illumination, too, but I’ll show you that. For now, pretend the illuminant is equal energy, too.)

The point is that in the real world, with real spectra — from pine needles or from some Munsell color chip or anything — every “metameric black” we compute could be added to any other spectrum we have, without affecting the color we perceive.

This, if northing else, reminds us that color is our perception of a light spectrum. And there are physical spectra in the visble region that appear black to us.

We can confirm that the fundamental f has precisely the same “R,G,B” values as the equal energy signal, by computing A’ f

Recall, or recompute A’ ee

The same. We have indeed split ee into two parts

ee = f + n,

one of which (n) contributes nothing to the final R,G,B values, the other of which (f) contributes everything. This is why f is called the fundamental.

And, our observation that the columns of E are projected onto themselves says that the columns of E are also fundamental spectra.

And they are a basis, so every fundamental — everything we see — could be written as a linear combination of the columns of E. They’re not an orthonormal basis, so we would have to use the dual basis to compute coefficients of the linear combinations. Well, that’s exactly what we accomplish by applying A’ (the dual basis) to the fundamental.

Those 3 numbers (A’f = A’ ee) we get describe the fundamental f in terms of the basis vectors E1, E2, E3:

f = R E1 + G E2 + B E3.

They are simply the components of f wrt the basis E1, E2, E3.

They are not, however, the coefficients of the original signal ee. That lies not entirely in the space spanned by E — that’s what the metameric black n is, the part of the original signal that is not in the space spanned by E.

(They also are not necessarily the RGB values we would use in RGB color space. After all, the components of f could be any real numbers, while RGB color space has restricted values, whether 0 to 1, 0 to 100, or 0 to 255. Don’t worry, we’ll get there.)

Summary

Cohen has a matrix A of “color matching functions”; it has 3 columns and lots more rows.

We construct a dual basis E via

.

Since “dual basis” is a dual relationship, we could have — and did! — compute A from a given E:

We construct two dispersion matrices

Me = E’E,

Ma = A’A,

and we saw that they were inverses.

We did Cholesky decompositions (LU decompositions of positive definite symmetric matrices) of Ma and Me, getting

Ga Ga’ = Ma,

Ge Ge’ = Me.

Then we used Ge and Ga to construct two orthonormal bases F1 and F2

F1 = A Ge,

F2 = E Ga.

Finally, Cohen computed the projection operator R onto the non-nullspace part of the domain of A’ (the preimage of the range of A’, equivalently the image of A).

Most importantly, at the end, we saw that an equal-energy (constant) spectrum could be split into two parts (f and n), one of which (f) contained all the “R,G,B” nformation, the other of which (n) was metameric black, invisible.

Next, I will show what I would do differently.

Color: HSB (HSV, HSI) again

I need to learn to write better summaries.

Here are the summaries I should have written for the HSB post, and for the tint-tone-shade post.

First, I am going to describe what “saturation” and “brightness- intensity- value” mean in the HSB (HSI, HSV) color space.

I cannot over emphasize that I am discussing what these words mean in a particular color space. As far as I know, these are not universal truths, but particular approaches to them. (I could be wrong: these definitions may be more general than I know, but that’s the point – I don’t know that they are universal.)

I still don’t like the non-technical definitions which I see all over the place, but I will admit that there is a place for them. Color is a physiological response, so many questions of color boil down to describing our response to it. And that is, ultimately, personal and subjective; we don’t all have the same eyes.

Someday I may find other quantitative definitions of “saturation” and “intensity”. This is what I have, so let me summarize it. (Of the three equivalent terms — brightness, intensity, value — I am going to use “intensity” in this post.)

Let me say that again with different emphasis. Brightness, intensity, and value are equivalent terms in these color spaces.

If you need or want far more than this summary, you could start by reading my original color posts.

I personally understand HSB by its relationship to RGB. Although hue, H, is the most complicated to describe mathematically, is the one that is obvious to me. All we do is choose an angle; keep it constant, and we’re working with tints, tones, and shades of a fixed color.

It is convenient, but not essential, to normalize R,G,B to the range 0..1.
If we’re given values in the range 0 to 255, just divide by 255. The equations work perfectly well if the range is 0 to 255; I just want to be able to refer to a numerical maximum, and I am choosing 1 instead of 255 because that is what my code used. We could just as well write the values as percentages, ranging from 0 to 100.

Intensity is computed easily.

• find the maximum value of R, G, B.
• that’s the intensity.

(It doesn’t matter if more than one of R,G,B has the maximum value; after all, it’s the same maximum. Intensity does not depend on which color has the maximum.)

We say that a color is fully intense if its intensity = 1; this is equivalent to: at least one of R,G,B = 1.

Saturation is computed as follows.

• Find the maximum and minimum values of R,G,B,
• subtract them: range = maximum – minimum
• and divide the range by the maximum.
• i.e. saturation = (maximum – minimum)/maximum.

It is worth noting that we can write the saturation as

saturation = 1 – minimum/maximum.

Among other things, that tells us that we can hold the saturation constant by holding constant the ratio of minimum to maximum.

There are interesting implications at the boundaries.

• Saturation = 1 if and only if minimum = 0.
• Saturation = 0 if and only if minimum = maximum.

These, in turn, translate to

• Saturation = 1 if and only if at least one of R,G,B = 0.
• Saturation = 0 iff R = G = B, i.e. we have gray (including possibly white or black).

(And that last property is why people sometimes say that saturation is “how far from gray” a color is.)

We say that a color is fully saturated if its saturation = 1; this is equivalent to one of R,G,B = 0.

So a color is fully saturated and fully intense if and only if one of R,G,B is 0 and one of R,G,B is 1.

(Incidentally, the HSB (HSI, HSV) definition of intensity does not capture the fact that a fully saturated and fully intense yellow appears more luminous than a fully intense and fully saturated blue.)

Second, let me summarize the relationship of intensity and saturation to tint-tone-shade.

I personally understand monochromatic variations of a fixed hue as tint-tone-shade. This does not help me mix colors, but it gives me a firm handle on variations of one color.

If we start with a fully intense and fully saturated hue, we reduce the saturation by increasing the minimum value of R,G,B; and this has no effect on intensity. Conceptually, that should be the equivalent of mixing it with white. And that’s what I showed in the earlier post: that reducing the saturation produces a tint.

Similarly, if we start with a fully intense and fully saturated hue, we reduce the intensity by decreasing the maximum value of R,G,B. This should be the equivalent of mixing it with black light. I know, that isn’t really possible, but it’s what’s happening: reducing the intensity produces a shade.

(If we start with a color which is not fully saturated, to reduce the intensity while holding saturation constant requires that we change both the maximum and minimum values of R,G,B. All we have to do is keep the ratio minimum/maximum constant.)

Finally, if we start with a fully intense and fully saturated hue, we produce a tone by decreasing the saturation (i.e. by increasing the minimum value of R,G,B) and by decreasing the intensity (i.e. by decreasing the maximum value of R,G,B).

Let me summarize. In HSB, HSV, or HSI:

• a tint is fully intense but not fully saturated;
• a shade is fully saturated but not fully intense;
• a tone is neither fully saturated nor fully intense.

I encourage you to play with the color picker on your computer. Let it show you the HSB and RGB values for tints, tones, and shades.