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:

g1102a

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:

g1102b

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 A^T\ — the (transpose of) the color matching functions A — to the reflected spectrum S, getting:

A^T\ S = \{356.817,\ 354.641,\ 271.109\}

I claim that these numbers are components of the fundamental of the reflected spectrum wrt the basis dual to A^T\ . 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 A^T\ — is the total cone response — the photopic response — of the standard observer.

Here’s the second row of A^T\ 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.

g1102c

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)…

T31 = \left(\begin{array}{ccc} 2.76888 & 1.75175 & 1.13016 \\ 1 & 4.5907 & 0.0601 \\ 0 & 0.05651 & 5.59427\end{array}\right)

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

XYZ = T31 RGB,

so

RGB = T31^{-1}\ XYZ = \{0.133246,\ 0.115616,\ 0.0903103\}\ .

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

E^T A = I\ .

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 E^T A\ .

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 A^T\ 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 A^T\ .

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 E^T dual to A — but that is the same as finding a basis E dual to A^T\ . The two equations

E^T A = I

and

A^T E = I

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

E^T A = I

E^T = A^{-1}\ .

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 (A^T A\ ), however, is invertible, and we construct the pseudo-inverse

(A'A)^{-1} A^T

and substitute it for A^{-1}\ , getting

E^T = A (A'A)^{-1}\ .

(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 E^T A\ . We might as well have found E^T first.

When we check it, we see that E^T A = I\ .

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

E = \left(\begin{array}{ccc} 0.000358806 & -0.000360371 & 0.000884809 \\ 0.00330737 & -0.00334373 & 0.00929769 \\ 0.0302172 & -0.0308581 & 0.0885709 \\ 0.0643748 & -0.0678399 & 0.241909 \\ 0.0116003 & -0.0190679 & 0.237446 \\ -0.0843071 & 0.0839105 & 0.126398 \\ -0.190772 & 0.216831 & 0.060478 \\ -0.329572 & 0.420346 & 0.0457603 \\ -0.286455 & 0.453811 & 0.0276004 \\ -0.0869933 & 0.321465 & -0.00307957 \\ 0.209686 & 0.0758441 & -0.0425853 \\ 0.438522 & -0.150742 & -0.0706912 \\ 0.418257 & -0.20126 & -0.0640185 \\ 0.232068 & -0.121142 & -0.0349616 \\ 0.0872125 & -0.0467666 & -0.0130635 \\ 0.0249132 & -0.0134701 & -0.00372503 \\ 0.00608984 & -0.00330711 & -0.000909677 \\ 0.00157141 & -0.000868448 & -0.000233816 \\ 0.000348975 & -0.000172577 & -0.0000531557 \\ 0.0000923158 & -0.0000402793 & -0.0000143874 \\ 0. & 0. & 0.\end{array}\right)

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

g1102d

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):

g1102e

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):

g1102f

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 A^T\ . And, since A^T\ 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 A^T\ 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 A^T\ 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 A^T\ 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 A^T\ 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 A^T\ , then a projection operator R onto that pre-image, then get the fundamental by applying the projection operator to the reflected spectrum.

First, let

A^T\ = u\ w\ v^T\ .

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

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

Here are the leftmost 3 columns of v:

v1 = \left(\begin{array}{ccc} -0.00201177 & 0.00143737 & 0.000374614 \\ -0.0210175 & 0.0149756 & 0.00335666 \\ -0.199597 & 0.142652 & 0.0305073 \\ -0.539263 & 0.385948 & 0.0611556 \\ -0.513526 & 0.367224 & -0.00344907 \\ -0.258102 & 0.159433 & -0.11971 \\ -0.121682 & -0.0133979 & -0.265642 \\ -0.134798 & -0.177324 & -0.481195 \\ -0.193185 & -0.30181 & -0.471062 \\ -0.246728 & -0.377774 & -0.254046 \\ -0.283178 & -0.410206 & 0.0979204 \\ -0.273373 & -0.377298 & 0.391838 \\ -0.201477 & -0.270573 & 0.408489 \\ -0.102016 & -0.135415 & 0.232377 \\ -0.0370629 & -0.0489598 & 0.0880787 \\ -0.0104736 & -0.0138137 & 0.0252276 \\ -0.00254533 & -0.00335417 & 0.00617542 \\ -0.000641283 & -0.000842035 & 0.00160261 \\ -0.000163267 & -0.00021855 & 0.000343642 \\ -0.0000487128 & -0.0000661713 & 0.0000876566 \\ 0. & 0. & 0.\end{array}\right)

We should look at those. First alone:

g1102g

… and then juxtaposed with the dual basis:
g1102g2

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 v1^T v1\ ; 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.

g1102h

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 A^T\ onto the preimage; we just take the product of v1 and v1^T in the other order:

R = v1\ v1^T

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

g1102 R

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 A^T\ by computing A^T\ n\ :

A^T\ n = \{-1.01618*10^-13,\ -1.09023*10^-13,\ 4.93919*10^-15\}

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

A^T\ S = \{356.817,\ 354.641,\ 271.109\}

for A^T\ applied to the reflected spectrum. If we now apply A^T\ to the fundamental, we get

A^T\ f = \{356.817,\ 354.641,\ 271.109\}

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.

g1102i

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.

g1102j

They are the same.

Thus we confirm that the 3 values A^T S 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.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: