G3dSpectrumImpl.mesa
Copyright Ó 1984, 1992 by Xerox Corporation. All rights reserved.
Glassner, May 12, 1989 9:45:55 am PDT
Bloomenthal, July 14, 1992 2:09 pm PDT
DIRECTORY G2dMatrix, G3dBasic, G3dSpectrum, IO, RealFns;
G3dSpectrumImpl: CEDAR PROGRAM
IMPORTS G2dMatrix, RealFns
EXPORTS G3dSpectrum
~ BEGIN
OPEN G3dSpectrum;
Type Declarations
SpectralSample:   TYPE ~ G3dSpectrum.SpectralSample;
Spectrum:     TYPE ~ G3dSpectrum.Spectrum;
SpectrumRep:   TYPE ~ G3dSpectrum.SpectrumRep;
RealSequence:   TYPE ~ G3dSpectrum.RealSequence;
RealSequenceRep:  TYPE ~ G3dSpectrum.RealSequenceRep;
Matrix:     TYPE ~ G3dSpectrum.Matrix;
RGB:      TYPE ~ G3dSpectrum.RGB;
XYZ:      TYPE ~ G3dSpectrum.XYZ;
Triple:     TYPE ~ G3dSpectrum.Triple;
Globals
xSequence, ySequence, zSequence:   RealSequence;
conSequence, sinSequence, cosSequence:  RealSequence;
conRGB, sinRGB, cosRGB:      RGB;
xBar, yBar, zBar:        Spectrum;
conSpectrum, sinSpectrum, cosSpectrum:  Spectrum;
weightsMat:          Matrix;
currentTransformMatrix:       Matrix;
ntscXYZtoRGBMatrix:       PUBLIC Matrix ¬ [
[ 1.967, -0.955, 0.065],
[-0.548, 1.938, -0.130],
[-0.297, -0.027, 0.982]
];
Basic Operations on a Single Spectrum
New: PUBLIC PROC [firstWavelength: REAL ¬ 380.0, deltaWavelength: REAL ¬ 30.0, numSamples: INT ¬ 16, amplitudes: RealSequence ¬ NIL] RETURNS [new: Spectrum] ~ {
new ¬ NEW[SpectrumRep[numSamples]];
new.length ¬ numSamples;
FOR i: CARDINAL IN [0 .. numSamples) DO
new[i].wavelength ¬ firstWavelength + i * deltaWavelength;
IF amplitudes = NIL
THEN new[i].amplitude ¬ 0.0
ELSE new[i].amplitude ¬ amplitudes[MIN[i, amplitudes.length-1]];
ENDLOOP;
};
Flatten: PUBLIC PROC [s: Spectrum, amplitude: REAL] RETURNS [Spectrum] ~ {
FOR i: INT IN [0 .. s.length) DO
s[i].amplitude ¬ amplitude;
ENDLOOP;
RETURN[s];
};
Normalize: PUBLIC PROC [s: Spectrum, totalEnergy: REAL ¬ 1.0] RETURNS [Spectrum] ~ {
energy: REAL ¬ 0.0;
FOR i: INT IN [0 .. s.length) DO energy ¬ energy + s[i].amplitude; ENDLOOP;
IF energy = 0.0 THEN RETURN [s];
energy ¬ totalEnergy / energy;
FOR i: INT IN [0 .. s.length) DO s[i].amplitude ¬ s[i].amplitude * energy; ENDLOOP;
RETURN[s];
};
Integrate: PUBLIC PROC [s: Spectrum] RETURNS [energy: REAL] ~ {
energy ¬ 0.0;
FOR i: INT IN [0 .. s.length) DO energy ¬ energy + s[i].amplitude; ENDLOOP;
RETURN[energy];
};
Scale: PUBLIC PROC [s: Spectrum, scale: REAL] RETURNS [Spectrum] ~ {
FOR i: INT IN [0 .. s.length) DO s[i].amplitude ¬ s[i].amplitude * scale; ENDLOOP;
RETURN[s];
};
Offset: PUBLIC PROC [s: Spectrum, offset: REAL] RETURNS [Spectrum] ~ {
FOR i: INT IN [0 .. s.length) DO s[i].amplitude ¬ s[i].amplitude + offset; ENDLOOP;
RETURN[s];
};
Duplicate: PUBLIC PROC [s: Spectrum] RETURNS [new: Spectrum] ~ {
new ¬ NEW[SpectrumRep[s.length]];
new.length ¬ s.length;
FOR i: INT IN [0 .. s.length) DO new[i] ¬ s[i]; ENDLOOP;
};
Sample: PUBLIC PROC [s: Spectrum, wavelength: REAL] RETURNS [REAL] ~ {
i, j: INT;
dist: REAL;
IF s.length = 0 THEN RETURN [0.0];
IF wavelength <= s[0].wavelength THEN RETURN [s[0].amplitude];
IF wavelength >= s[s.length-1].wavelength THEN RETURN [s[s.length-1].amplitude];
i ¬ 0; WHILE s[i].wavelength < wavelength DO i ¬ i+1; ENDLOOP;
j ¬ i-1;
dist ¬ s[i].wavelength - s[j].wavelength;
IF dist = 0.0 THEN {
j: INT ¬ i-1;
WHILE j>0 AND s[j].wavelength = s[i].wavelength DO j ¬ j-1; ENDLOOP;
IF j=0 AND s[j].wavelength = s[i].wavelength THEN RETURN[s[0].amplitude];
};
dist ¬ (wavelength - s[j].wavelength)/dist;
RETURN [s[j].amplitude + dist *(s[i].amplitude - s[j].amplitude)];
};
Color Space Conversions
SetTransformMatrix: PUBLIC PROC [transform: Matrix ¬ ntscXYZtoRGBMatrix] ~ {
currentTransformMatrix ¬ transform;
conRGB ¬ SpectrumToRGB[conSpectrum];
sinRGB ¬ SpectrumToRGB[sinSpectrum];
cosRGB ¬ SpectrumToRGB[cosSpectrum];
weightsMat ¬ [
[conRGB.R, conRGB.G, conRGB.B],
[sinRGB.R, sinRGB.G, sinRGB.B],
[cosRGB.R, cosRGB.G, cosRGB.B]
];
weightsMat ¬ G2dMatrix.Invert[weightsMat];
};
TransformFromChromaticities: PUBLIC PROC [r, g, b, w: Triple] RETURNS [transform: Matrix] ~ {
k: Matrix ¬ [r, g, b];
wnew: Triple ¬ IF w.y # 0.0 THEN [w.x/w.y, 1.0, w.z/w.y] ELSE [0.0, 1.0, 0.0];
ki: Matrix ¬ G2dMatrix.Invert[k];
v: Triple ¬ Pt3Mul[wnew, ki];
gm: Matrix ¬ [[v.x, 0.0, 0.0], [0.0, v.y, 0.0], [0.0, 0.0, v.z]];
n: Matrix ¬ G2dMatrix.Mul[gm, k];
transform ¬ G2dMatrix.Invert[n];
};
SpectrumToRGB: PUBLIC PROC [s: Spectrum] RETURNS [rgb: RGB] ~ {
xyz: XYZ ¬ [0.0, 0.0, 0.0];
news: Spectrum ¬ Conform[s, xBar];
t: Triple;
FOR i: INT IN [0 .. xBar.length) DO
xyz.x ¬ xyz.x + (s[i].amplitude * xBar[i].amplitude);
xyz.y ¬ xyz.y + (s[i].amplitude * yBar[i].amplitude);
xyz.z ¬ xyz.z + (s[i].amplitude * zBar[i].amplitude);
ENDLOOP;
t ¬ Pt3Mul[[xyz.x, xyz.y, xyz.z], currentTransformMatrix];
rgb ¬ [t.x, t.y, t.z];
};
RGBToSpectrum: PUBLIC PROC [rgb: RGB] RETURNS [Spectrum] ~ {
t: Triple ¬ Pt3Mul[[rgb.R, rgb.G, rgb.B], weightsMat];
RETURN[Blend[cosSpectrum, Blend[conSpectrum, sinSpectrum, t.x, t.y], t.z, 1.0]];
};
Pt3Mul: PROC [v: Triple, m: Matrix] RETURNS [Triple] ~ {
RETURN[[
(v.x * m.row1.x) + (v.y * m.row2.x) + (v.z * m.row3.x),
(v.x * m.row1.y) + (v.y * m.row2.y) + (v.z * m.row3.y),
(v.x * m.row1.z) + (v.y * m.row2.z) + (v.z * m.row3.z)
]];
};
Basic Operations on Two Spectra
Union: PUBLIC PROC [a, b: Spectrum] RETURNS [newa, newb: Spectrum] ~ {
swap, ulen: INT ¬ a.length;
union: RealSequence ¬ NEW[RealSequenceRep[a.length + b.length]];
first build the union of the two wavelength lists
FOR i: INT IN [0 .. a.length) DO union[i] ¬ a[i].wavelength; ENDLOOP;
FOR i: INT IN [0 .. b.length) DO
unused: BOOLEAN ¬ TRUE;
FOR j: INT IN [0 .. a.length) DO
IF a[j].wavelength = b[i].wavelength THEN unused ¬ FALSE;
ENDLOOP;
IF unused THEN {
union[ulen] ¬ b[i].wavelength;
ulen ¬ ulen + 1;
};
ENDLOOP;
now bubble-sort the wavelength list (not a bad idea for such a small list)
FOR i: INT IN [0 .. ulen) DO
swap ¬ -1;
FOR j: INT IN [i+1 .. ulen) DO
IF union[j] < union[i] THEN swap ¬ j;
ENDLOOP;
WHILE swap >= 0 DO
tmp: REAL;
tmp ¬ union[i]; union[i] ¬ union[swap]; union[swap] ¬ tmp;
swap ¬ -1;
FOR j: INT IN [i+1 .. ulen) DO
IF union[j] < union[i] THEN swap ¬ j;
ENDLOOP;
ENDLOOP;
ENDLOOP;
and finally build the composite spectra
newa ¬ NEW[SpectrumRep[ulen]];
newb ¬ NEW[SpectrumRep[ulen]];
newa.length ¬ newb.length ¬ ulen;
FOR i: INT IN [0 .. ulen) DO
newa[i].wavelength ¬ newb[i].wavelength ¬ union[i];
newa[i].amplitude ¬ Sample[a, union[i]];
newb[i].amplitude ¬ Sample[b, union[i]];
ENDLOOP;
};
Conform: PUBLIC PROC [a, b: Spectrum] RETURNS [newa: Spectrum] ~ {
newa ¬ NEW[SpectrumRep[b.length]];
newa.length ¬ b.length;
FOR i: INT IN [0 .. b.length) DO
newa[i].wavelength ¬ b[i].wavelength;
newa[i].amplitude ¬ Sample[a, b[i].wavelength];
ENDLOOP;
};
Add: PUBLIC PROC [a, b: Spectrum] RETURNS [c: Spectrum] ~ {
newa, newb: Spectrum;
[newa, newb] ¬ Union[a, b];
c ¬ NEW[SpectrumRep[newa.length]];
c.length ¬newa.length;
FOR i: INT IN [0 .. a.length) DO
c[i].wavelength ¬ newa[i].wavelength;
c[i].amplitude ¬ newa[i].amplitude + newb[i].amplitude;
ENDLOOP;
};
Sub: PUBLIC PROC [a, b: Spectrum] RETURNS [c: Spectrum] ~ {
newa, newb: Spectrum;
[newa, newb] ¬ Union[a, b];
c ¬ NEW[SpectrumRep[newa.length]];
c.length ¬newa.length;
FOR i: INT IN [0 .. a.length) DO
c[i].wavelength ¬ newa[i].wavelength;
c[i].amplitude ¬ newa[i].amplitude - newb[i].amplitude;
ENDLOOP;
};
Lerp: PUBLIC PROC [a, b: Spectrum, alpha: REAL] RETURNS [c: Spectrum] ~ {
newa, newb: Spectrum;
[newa, newb] ¬ Union[a, b];
c ¬ NEW[SpectrumRep[newa.length]];
c.length ¬newa.length;
FOR i: INT IN [0 .. a.length) DO
c[i].wavelength ¬ newa[i].wavelength;
c[i].amplitude ¬ newa[i].amplitude + alpha * (newb[i].amplitude - newa[i].amplitude);
ENDLOOP;
};
Blend: PUBLIC PROC [a, b: Spectrum, aScl, bScl: REAL] RETURNS [c: Spectrum] ~ {
newa, newb: Spectrum;
[newa, newb] ¬ Union[a, b];
c ¬ NEW[SpectrumRep[newa.length]];
c.length ¬newa.length;
FOR i: INT IN [0 .. a.length) DO
c[i].wavelength ¬ newa[i].wavelength;
c[i].amplitude ¬ (aScl * newa[i].amplitude) + (bScl * newb[i].amplitude);
ENDLOOP;
};
Filter: PUBLIC PROC [a, b: Spectrum] RETURNS [c: Spectrum] ~ {
newa, newb: Spectrum;
[newa, newb] ¬ Union[a, b];
c ¬ NEW[SpectrumRep[newa.length]];
c.length ¬newa.length;
FOR i: INT IN [0 .. a.length) DO
c[i].wavelength ¬ newa[i].wavelength;
c[i].amplitude ¬ newa[i].amplitude * newb[i].amplitude;
ENDLOOP;
};
Damp: PUBLIC PROC [a, b: Spectrum] RETURNS [c: Spectrum] ~ {
newa, newb: Spectrum;
[newa, newb] ¬ Union[a, b];
c ¬ NEW[SpectrumRep[newa.length]];
c.length ¬newa.length;
FOR i: INT IN [0 .. a.length) DO
c[i].wavelength ¬ newa[i].wavelength;
IF newb[i].amplitude # 0.0
THEN c[i].amplitude ¬ newa[i].amplitude / newb[i].amplitude
ELSE c[i].amplitude ¬ 0.0;
ENDLOOP;
};
Start Code
CIE 1931 2-degree observer matching functions from 380 to 780 nm in 5nm steps
xSamples: ARRAY [0 .. 80] OF REAL ¬ [
0.0014, 0.0022, 0.0042, 0.0076, 0.0143, 0.0232, 0.0435, 0.0776, 0.1344,
0.2148, 0.2839, 0.3285, 0.3483, 0.3481, 0.3362, 0.3187, 0.2908, 0.2511,
0.1954, 0.1421, 0.0956, 0.0580, 0.0320, 0.0147, 0.0049, 0.0024, 0.0093,
0.0291, 0.0633, 0.1096, 0.1655, 0.2257, 0.2904, 0.3597, 0.4334, 0.5121,
0.5945, 0.6784, 0.7621, 0.8425, 0.9163, 0.9786, 1.0263, 1.0567, 1.0622,
1.0456, 1.0026, 0.9384, 0.8544, 0.7514, 0.6424, 0.5419, 0.4479, 0.3608,
0.2835, 0.2187, 0.1649, 0.1212, 0.0874, 0.0636, 0.0468, 0.0329, 0.0227,
0.0158, 0.0114, 0.0081, 0.0058, 0.0041, 0.0029, 0.0020, 0.0014, 0.0010,
0.0007, 0.0005, 0.0003, 0.0002, 0.0002, 0.0001, 0.0001, 0.0000, 0.0000
];
ySamples: ARRAY [0 .. 80] OF REAL ¬ [
0.0000, 0.0001, 0.0001, 0.0002, 0.0004, 0.0006, 0.0012, 0.0022, 0.0040,
0.0073, 0.0116, 0.0168, 0.0230, 0.0298, 0.0380, 0.0480, 0.0600, 0.0739,
0.0910, 0.1126, 0.1390, 0.1693, 0.2080, 0.2586, 0.3230, 0.4073, 0.5030,
0.6082, 0.7100, 0.7932, 0.8620, 0.9149, 0.9540, 0.9803, 0.9950, 1.0002,
0.9950, 0.9786, 0.9520, 0.9154, 0.8700, 0.8163, 0.7570, 0.6949, 0.6310,
0.5668, 0.5030, 0.4412, 0.3810, 0.3210, 0.2650, 0.2170, 0.1750, 0.1382,
0.1070, 0.0816, 0.0610, 0.0446, 0.0320, 0.0232, 0.0170, 0.0119, 0.0082,
0.0057, 0.0041, 0.0029, 0.0021, 0.0015, 0.0010, 0.0007, 0.0005, 0.0004,
0.0003, 0.0002, 0.0001, 0.0001, 0.0001, 0.0000, 0.0000, 0.0000, 0.0000
];
zSamples: ARRAY [0 .. 80] OF REAL ¬ [
0.0065, 0.0105, 0.0201, 0.0362, 0.0679, 0.1102, 0.2074, 0.3713, 0.6456,
1.0391, 1.3856, 1.6230, 1.7471, 1.7826, 1.7721, 1.7441, 1.6692, 1.5281,
1.2876, 1.0419, 0.8130, 0.6162, 0.4652, 0.3533, 0.2720, 0.2123, 0.1582,
0.1117, 0.0782, 0.0573, 0.0422, 0.0298, 0.0203, 0.0134, 0.0087, 0.0057,
0.0039, 0.0027, 0.0021, 0.0018, 0.0017, 0.0014, 0.0011, 0.0010, 0.0008,
0.0006, 0.0003, 0.0002, 0.0002, 0.0001, 0.0000, 0.0000, 0.0000, 0.0000,
0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000
];
xSequence ¬ NEW[RealSequenceRep[81]];
ySequence ¬ NEW[RealSequenceRep[81]];
zSequence ¬ NEW[RealSequenceRep[81]];
FOR i: INT IN [0 .. 81) DO
xSequence[i] ¬ xSamples[i];
ySequence[i] ¬ ySamples[i];
zSequence[i] ¬ zSamples[i];
ENDLOOP;
xBar ¬ New[firstWavelength:380, deltaWavelength: 5, numSamples: 81, amplitudes: xSequence];
yBar ¬ New[firstWavelength:380, deltaWavelength: 5, numSamples: 81, amplitudes: ySequence];
zBar ¬ New[firstWavelength:380, deltaWavelength: 5, numSamples: 81, amplitudes: zSequence];
now set up for RGB->Spectrum conversions with the three necessary spectra
currentTransformMatrix ¬ ntscXYZtoRGBMatrix;
conSequence ¬ NEW[RealSequenceRep[81]];
sinSequence ¬ NEW[RealSequenceRep[81]];
cosSequence ¬ NEW[RealSequenceRep[81]];
FOR i: INT IN [0 .. 81) DO
theta: REAL ¬ 6.283185*i/80.0;
conSequence[i] ¬ 1.0;
sinSequence[i] ¬ RealFns.Sin[theta];
cosSequence[i] ¬ RealFns.Cos[theta];
ENDLOOP;
conSpectrum ¬ New[firstWavelength:380, deltaWavelength: 5, numSamples: 81, amplitudes: conSequence];
sinSpectrum ¬ New[firstWavelength:380, deltaWavelength: 5, numSamples: 81, amplitudes: sinSequence];
cosSpectrum ¬ New[firstWavelength:380, deltaWavelength: 5, numSamples: 81, amplitudes: cosSequence];
conRGB ¬ SpectrumToRGB[conSpectrum];
sinRGB ¬ SpectrumToRGB[sinSpectrum];
cosRGB ¬ SpectrumToRGB[cosSpectrum];
weightsMat ¬ [
[conRGB.R, conRGB.G, conRGB.B],
[sinRGB.R, sinRGB.G, sinRGB.B],
[cosRGB.R, cosRGB.G, cosRGB.B]
];
weightsMat ¬ G2dMatrix.Invert[weightsMat];
END.