<<>>
<<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.

<<>>