|
1
|
On this web site, the
phrase "vectorial color" denotes a set of methods in which lights are
represented by
vectors in the color space developed by Jozef Cohen. The vectors are
tristimulus vectors, based on
the orthonormal color matching functions. The goal is to make better use of existing
color-matching data. Whether the 10° data are better
than the 2° data, or new
experimental data are needed, those are separate questions. You choose
your data,
and the vectorial methods will help to make use of them.
The starting point for
vectorial color is to take existing color
matching functions, such as
the 2° observer
functions, and add and subtract them to make orthonormal color matching
functions. The "adding and subtracting" is done by matrix algebra and
in one step uses the well-known Gram-Schmidt method of
orthonormalization. The purpose of this web page (as of 2006 June) is
to present a few code fragments and short working routines, to show how
the work is done. With regard to some methods, such as Gram-Schmidt, it
is easier to write a program than to describe the method in words and
symbols.
I do all my calculations in a science-programming language called O-Matrix, which
allows matrix operations to be written as easily as operations on
scalars. If a and b are matrices previously created,
one can write
c = a*b
to right-multiply a by b. The matrix c need not be declared or
dimensioned in
advance. In short, steps that are simple in concept remain simple in
the
computer program. Another language similar to O-Matrix is Matlab. The
O-Matrix manual exists only as web pages, and you can view a
version on their
web site.
If you would like to discuss color calculations, or request more
details, please contact jim@jimworthey.com
, or call 301-977-3551 in the USA.
|
|
|
|
2
|
Suppose for a moment that
the orthonormal color matching functions have been found. For the 2° observer, they
are as graphed above. Conceptually, they can be 3 column matrices,
ω1, ω2, ω3. The three column vectors can be grouped together in an N×3 matrix, where N is the
number of wavelengths:
Ω = [ω1 ω2 ω3] .
This
is the kind of notation used in the articles, with the Greek letters,
big and little omega.
For programming, it is convenient to call Ω something like
OrthoBasis:
Ω →
OrthoBasis .
If
an individual vector is needed, it can be referred to in O-Matrix as
OrthoBasis.col(j), where j = 1, 2, or 3 as needed. The
notation ".col(2)" selects the second column of a matrix, for example.
Now suppose that D65 is a column vector representing the spectrum of
standard daylight, D65. We wish to calculate its tristimulus vector,
which can be called V65. Then
V65 =
OrthoBasis'*D65 .
That
line of code is in a larger font to be sure that you can see the
apostrophe after OrthoBasis. It indicates matrix transpose. Using the
proper notation then, the tristimulus vector of any light can be found
in one line of computer code. In order to do the operations of
colorimetry by matrix operations, one must align the functions of
wavelength properly: they must have the same number of elements and the
same wavelength domain.
|
 |
|
|
3
|
In
the visuals of the 2004
talk, "Color
Matching with Amplitude not Left Out," or in the new unpublished
article, "Vectorial Color," the idea of orthonormal color matching
functions is developed step-by-step. It is logical to have one
all-positive function, and to let
that be an achromatic function that is proportional to the usual y-bar.
But the achromatic function can be normalized. Suppose that ybar is
usual 2° observer
function. Then a code fragment will show the steps:
prod = ybar'*ybar
omega1 =
ybar/sqrt(prod) .
Again the apostrophe
denotes matrix transpose. sqrt(prod) is then the vector length of ybar.
Dividing ybar by its length gives omega1, a column matrix proportional
to ybar,
but with length = 1. To verify,
print
omega1'*omega1 .
The
result should
equal one.
Now suppose that red is
a
red cone sensitivity function. We know that omega1 is a sum of red and
green, with certain proportions. To make an opponent function
orthnormal to omega1, we subtract from red the projection of omega1
onto
red, then normalize:
omega2 = red -
(red'*omega1)*omega1
prod =
omega2'*omega2
omega2 =
omega2/sqrt(prod)
In normal programming
fashion, the array omega2 holds an intermediate matrix, then the final
answer. The steps above illustrate the method of Gram-Schmidt
orthonormalization. The following O-Matrix routine calculates an
orthonormal set from the columns of any matrix:
GramSchmidt(given, orthonorm)
.
The link will bring
up a text file containing the program. To run it, you must save the
text in a file called GramSchmidt.oms. Calling the Gram-Schmidt routine
by the following steps will create OrthoBasis:
given = [ybar,
red, blue]
NumCols = 3
NumRows =
rowdim(given)
OrthoBasis =
zeros(NumRows,NumCols)
GramSchmidt(given,
OrthoBasis)
It is assumed that
there is a set of cone functions {red, green, blue} that are linear
combinations of the usual xbar, ybar, zbar of the 2° observer, and
indeed that
ybar is a linear combination of red and green only. For background,
see the Vectorial Color
manuscript,
especially the appendices and Figure 1.
Also see item 4 below.
|
 |
|
|
4
|
Just above, and in Figure 1 of the Vectorial
Color Manuscript, reference is made to cone sensitivity functions
that are linear combinations of the usual 2° observer. I take
the idea from Guth, who borrowed it from Smith and Pokorny. Guth's
formula for the cone functions is recalled in Appendix A of the "Render Calc" article.
Taking the matrix
transpose of that version, including a transposed version of the square
matrix, we can write in O-Matrix:
M1 = { [0.2435,
-0.3954, 0], ...
[0.8524,
1.1642, 0], ...
[-0.0516, 0.0837,
0.6225] }
rgb = xyzbar*M1
The columns of xyzbar are
assumed to be the 2° observer
functions. That is xyzbar = [xbar, ybar, zbar] . The actual data are here, or can
be obtained from www.cvrl.org
. For my work, I store most spectral data in a particular format that
is convenient, but here the xyzbar data are lumped in one file,
just to show that they're available.
To get Guth's 1980 opponent model directly from xyzbar,
M3
= {[0, 0.7401, -0.0061], ...
[0.9341, -0.6801, -0.0212], ...
[0, -0.1567, 0.0314] }
atd = xyzbar*M3
OrthoBasis =
zeros(NumRows,NumCols)
GramSchmidt(atd,
OrthoBasis)
The
result is the same orthonormal basis as in item 3. The result should
also check against this
file.
|
 |
|
|
5
|
Matrix R.
Let A be a matrix whose columns are three linearly independent vectors,
such as three color matching functions. For example, A could be the
usual 2° observer
functions.
Then in O-Matrix, projector matrix R can be found by Cohen's formula:
A = xyzbar
R = A*inv(A'*A)*A'
In the special case that the columns
of A are an orthonormal set, for example if A = OrthoBasis, then A'*A
is an identity matrix. Its inverse is also an identity matrix and one
may write for this special case only,
R = OrthoBasis*OrthoBasis'
An O-Matrix function can
be used:
Matrix R is a large
matrix. For example, if A indeed contains the 2°
observer functions at 1 nm intervals and wavelength domain 360 to 830
nm, then R has dimensions 471×471, meaning that it
has 221841 elements and in double precision (8 bytes per floating point
number), it occupies 1774728 bytes. On modern PC hardware, this storage
space and the time to calculate R are absolutely not a concern.
However, it is not helpful to print out the matrix on paper. On my PC,
the time to calculate a 471×471 matrix R by
function RCohen() is about 13 ms. |
 |
|
|
6
|
Putting Matrix R to work. The
projector matrix R can be viewed in 2 ways. For routine use of the
orthonormal basis, matrix R is not needed at all. Matrix R can be
considered as background information, as the proof that the
locus of unit monochromats has an invariant shape, independent of the
coordinate system. The orthonormal basis is a shortcut for creating
color vectors without using matrix R itself.
On the other hand, matrix R can be a practical tool for curve-fitting
problems that arise in relation to color. If matrix A
contains any 1, 2, 3, or more linearly independent vectors, then the
projector matrix R can be found. Then if any vector X (with the same
number of rows as A or R) is multiplied by R, the result is the
projection of X into the column space of A. That is, Xstar = R*X, where
Xstar is the projection. An alternate wording is that Xstar is the
least-squares best fit to X, using the columns of A.
On the right are
graphed the red, green, and blue
sensitivities of a digital camera. In analyzing and applying this color
sensor, one might wish to find a linear combination of the camera's red
and green sensitivities that best approximates the human achromatic
sensitivity. An alternate statement of the same goal is that we wish to
project the human achromatic function into the the vector space of the
camera's red and green sensors. Let the camera's sensitivities be
called redcam, greencam, bluecam. Then
Arg = [redcam,
greencam]
Rrg = RCohen(Arg)
achcam =
Rrg*OrthoBasis.col(1)
Column matrix achcam is then the best approximation to the
human achromatic function, OrthoBasis.col(1), by a
combination of redcam and greencam. The
lower graph shows achcam compared to OrthoBasis.col(1).
More steps would be needed for the complete analysis of a camera. At
this moment, in 2006 June, I am preparing a paper on that topic, with a
co-author. For now, the topic is calculation. The methods on this page
can be used for such steps as these:
- Find an
orthonormal basis for the camera itself. By combining those functions,
find the
camera's locus of unit monochromats.
- Find the best fit to
all 3 human sensitivities by the camera sensitivities.
The actual
programming is reduced to a few steps, and the problem becomes one of
concepts: what do you want the camera to do?
|
|
|
|
7
|
Blackbody spectra can be
calculated by this routine:
blackbody(Tkelvin,
descrip, LamLoHi, vec) .
Tkelvin = Kelvin
temperature of the blackbody
descrip = a character
array that will be set with a description of this light
LamLoHi = the domain
of wavelength, such as [360, 830].
vec = the result, a
column matrix containing the blackbody spectrum
The blackbody
spectrum is adjusted so that it = 1 at 555 nm. In this case,
LamLoHi is an input, set by the calling program. My programs are
inconsistent about this, whether LamLoHi must be set by the calling
routine. Since I want all spectra to have the same domain, the domain
is often set by a global variable.
|
 |
|
|
|
|
|
|
|
|
|
|
|
