---

Spectrum.cc

Go to the documentation of this file.

/*
    File:       Spectrum.cc

    Function:   Provides classes for manipulating spectra

    Author:     Andrew Willmott

    Notes:      

*/

#include <stdlib.h>
#include "Spectrum.h"
#include "ColourSystem.h"
#include "VecUtil.h"

const Int   Spectrum::kComponents = 100;    // number of basis functions
const Float Spectrum::kStart = 350.0;       // start of spectrum
const Float Spectrum::kEnd = 800.0;         // end of spectrum

// CIE colour matching functions, x_bar, y_bar and z_bar
// 360-830nm in steps of 5nm
// good source for this stuff is http://cvision.ucsd.edu/basicindex.htm

const Float kCIE_Start = 360.0;
const Float kCIE_End = 830.0;
const Float kCIE_Delta = 5.0;

static Float tCIE_xbar[95] =
{
#include "cie_xbar.spectrum"    
};
RegularSpectralCurve kCIE_xbar(tCIE_xbar, 360, 830, 95);

static Float tCIE_ybar[95] =
{
#include "cie_ybar.spectrum"    
};
RegularSpectralCurve kCIE_ybar(tCIE_ybar, 360, 830, 95);

static Float tCIE_zbar[95] =
{
#include "cie_zbar.spectrum"    
};
RegularSpectralCurve kCIE_zbar(tCIE_zbar, 360, 830, 95);

Colour Spectrum::ToCIE_XYZ()
{
    Int     i;
    Colour  c = vl_0;
    
    for (i = 0; i < kComponents; i++)
    {
        ClrReal x = (i + 0.5) / kComponents;
        ClrReal lambda = mix(kStart, kEnd, x);
        
        c[0] += kCIE_xbar.Sample(lambda) * SELF[i];
        c[1] += kCIE_ybar.Sample(lambda) * SELF[i];
        c[2] += kCIE_zbar.Sample(lambda) * SELF[i];
    }

    return(c);
}

// --- Spectral curve methods -------------------------------------------------

SpectralCurve::~SpectralCurve()
{
}

SpectralCurve::operator Spectrum()
{
    Spectrum    result;
    Float       lambda, x;
    Int         i;

#ifdef USE_INT
    Float       delta, integral, lastIntegral;
    cout << "INTEGRATING" << endl;
    delta = (Spectrum::kEnd - Spectrum::kStart) / Spectrum::kComponents;
    lambda = Spectrum::kStart;
    lastIntegral = Integral(lambda);
    
    for (i = 0; i < Spectrum::kComponents; i++)
    {
        lambda += delta;
        integral = Integral(lambda);
        // normalised box functions, so we divide by the support size
        result[i] = (integral - lastIntegral) / delta;
        lastIntegral = integral;
    }
#else

    for (i = 0; i < Spectrum::kComponents; i++)
    {
        x = (i + 0.5) / Spectrum::kComponents;
        lambda = Spectrum::kStart + x * (Spectrum::kEnd - Spectrum::kStart);

        result[i] = Sample(lambda);
    }

#endif
    
    return(result);
}

Float SpectralCurve::Integral(Float lambda) const
{
    const Float kIntegrateDelta = 1.0;
    Float       s, sum = 0.0;
    Int         n = 0;
    
    for (s = Spectrum::kStart; s < lambda; s += kIntegrateDelta)
    {
        sum += Sample(s);
        n++;
    }
    
    if (n > 0)
        sum *= (lambda - Spectrum::kStart) / n;

    return(sum);
}

Float SpectralCurve::MaxCmpt() const
{
    const Float kMaxDelta = 1.0;
    Float       s, max = Sample(Spectrum::kStart);
    
    for (s = Spectrum::kStart + kMaxDelta; s <= Spectrum::kEnd; s += kMaxDelta)
    {
        max = Max(max, Sample(s));
    }
    
    return(max);
}

Void SpectralCurve::MakeImage(Image &image, ClrReal maxAmplitude)
// Plot spectrum into image
{
    Int     i, j;
    Float   *samples = new Float[image.Width()], ampMax = 0.0;
    Colour  *waveColour = new Colour[image.Width()];

    image.Clear(cBlack);

    for (i = 0; i < image.Width(); i++)
    {
        Float   x = (i + 0.5) / image.Width();
        Float   lambda = mix(Spectrum::kStart, Spectrum::kEnd, x);
        
        samples[i] = Sample(lambda);
        waveColour[i] = WavelengthToXYZ(lambda);
        waveColour[i] = csDisplay.XYZToGamutClip(waveColour[i]);
        ampMax = Max(samples[i], ampMax);
    }

    if (maxAmplitude < 0.0)
        maxAmplitude = ampMax;

    for (j = 0; j < image.Height(); j++)
    { 
        Float y = maxAmplitude * (j + 0.5) / image.Height();
        
        for (i = 0; i < image.Width(); i++)
            if (samples[i] >= y)
                image.SetPixel(i, j, waveColour[i]);
    }

    delete[] samples;
    delete[] waveColour;
}

Colour SpectralCurve::ToCIE_XYZ()
{
    Float   lambda, sample;
    Colour  c = vl_0;
    
    for (lambda = kCIE_Start; lambda <= kCIE_End; lambda += kCIE_Delta)
    {
        sample = Sample(lambda);
        c[0] += kCIE_xbar.Sample(lambda) * sample;
        c[1] += kCIE_ybar.Sample(lambda) * sample;
        c[2] += kCIE_zbar.Sample(lambda) * sample;
    }

    return(c);
}

static Float *MakeAccumData(Float *data, Int divisions)
{
    Int     i;
    Float   *result, sum = 0;

    result = new Float[divisions];

    for (i = 0; i < divisions; i++)
    {
        sum += data[i];
        result[i] = sum;
    }

    return(result);
}

RegularSpectralCurve::RegularSpectralCurve(Float *sdata, Int start, Int end, Int divisions) : 
    data(sdata),
    waveDivs(divisions),
    waveStart(start),
    waveEnd(end)
{
    waveDelta = (end - start) / (divisions - 1);
    accumData = MakeAccumData(data, divisions);
}

RegularSpectralCurve::~RegularSpectralCurve()
{
    delete[] accumData;
}

Float RegularSpectralCurve::Sample(Float lambda) const
// interpolate between two closest samples
{
    Float   result, x, dx;
    Int     bin1, bin2;

    if (lambda < waveStart || lambda > waveEnd)
        return(0.0);

    x = (lambda - waveStart) / waveDelta;
        
    bin1 = (Int) x;
    dx = x - bin1;
    if (dx > 1e-8)
        bin2 = bin1 + 1;
    else
        bin2 = bin1;

    CheckRange(bin1, 0, waveDivs, "Bad bin number");
    CheckRange(bin2, 0, waveDivs, "Bad bin number");

    result = mix(data[bin1], data[bin2], dx);

    return(result);
}

#ifdef UNFINISHED
Float RegularSpectralCurve::Integral(Float lambda) const
// interpolate between two closest samples
{
    Float   result, x, dx;
    Int     bin1, bin2;

    if (lambda < waveStart)
        return(0.0);
    else if (lambda > waveEnd)
        return(waveDelta * accumData[waveDivs - 1]);
    
    x = (lambda - waveStart) / waveDelta;
        
    bin1 = (Int) x;
    dx = x - bin1;
    if (dx > 1e-8)
        bin2 = bin1 + 1;
    else
        bin2 = bin1;

    CheckRange(bin1, 0, waveDivs, "Bad bin number");
    CheckRange(bin2, 0, waveDivs, "Bad bin number");

    result = waveDelta * mix(accumData[bin1], accumData[bin2], dx);

    return(result);
}
#endif

IrregularSpectralCurve::IrregularSpectralCurve(Float *amps, Float *wavs, Int divs) :
    amplitudes(amps),
    wavelengths(wavs),
    divisions(divs)
{
    accumAmplitudes = MakeAccumData(amps, divs);
}

IrregularSpectralCurve::IrregularSpectralCurve() :
    amplitudes(0),
    wavelengths(0),
    divisions(0),
    accumAmplitudes(0)
{
}

IrregularSpectralCurve::~IrregularSpectralCurve()
{
    delete[] accumAmplitudes;
}

IrregularSpectralCurve &IrregularSpectralCurve::ReadSDFile(FileName filename)
{
    String          line;
    StrConstArray   tokens;
    ifstream        s;
    Int             i = 0;
    
    cout << "parsing " << filename.GetPath() << endl;

    s.open(filename.GetPath());

    if (!s)
    {
        cerr << "(ReadSDFile) Cannot access " << filename.GetPath() << endl;
        return(SELF);
    }

    while (s)
    {
        if (line.ReadLine(s))
        {
            Split(line, tokens);
            if (tokens.NumItems() == 0 || tokens[0][0] == '#')
                ;
            else if (tokens.NumItems() == 1)
            {
                divisions = atoi(tokens[0]);
                cout << "setting size to " << divisions << endl;
                delete[] wavelengths;
                delete[] amplitudes;
                wavelengths = new Float[divisions];
                amplitudes = new Float[divisions];
                i = 0;
            }
            else if (tokens.NumItems() >= 2)
            {
                if (!wavelengths)
                    _Error("(ReadSDFile) didn't find number of divisions");
                wavelengths[i] = atof(tokens[0]);
                amplitudes[i] = atof(tokens[1]);
                i++;
            }
        }
    }   

    s.close();
    accumAmplitudes = MakeAccumData(amplitudes, divisions);
    
    return(SELF);
}

Float IrregularSpectralCurve::Sample(Float lambda) const
// interpolate between two closest samples
// need binary search or something to speed this up. (cache?)
// currently we extend the endpoint samples out to infinity either
// side. Perhaps should ramp down over X nm instead?
{
    Float   x, result;
    Int     i;

    if (divisions == 0)
        return(0.0);
    else if (divisions == 1)
        return(amplitudes[0]);

    if (lambda < wavelengths[0])
        return(amplitudes[0]);

    i = 1;
    while (i < divisions)
        if (lambda < wavelengths[i])
        {
            x = (lambda - wavelengths[i - 1])
                    / (wavelengths[i] - wavelengths[i - 1]);

            result = mix(amplitudes[i - 1], amplitudes[i], x);

            return(result);
        }
        else
            i++;

    return(amplitudes[divisions - 1]);
}

#ifdef UNFINISHED
Float IrregularSpectralCurve::Integral(Float lambda) const
// interpolate between two closest samples
// need binary search or something to speed this up. (cache?)
// currently we extend the endpoint samples out to infinity either
// side. Perhaps should ramp down over X nm instead?
{
    Float   x, result;
    Int     i;

    if (divisions == 0)
        return(0.0);
    else if (divisions == 1)
        return(amplitudes[0]);

    if (lambda < wavelengths[0])
        return(amplitudes[0]);

    i = 1;
    while (i < divisions)
        if (lambda < wavelengths[i])
        {
            x = (lambda - wavelengths[i - 1])
                    / (wavelengths[i] - wavelengths[i - 1]);

            result = mix(amplitudes[i - 1], amplitudes[i], x);

            return(result);
        }
        else
            i++;

    return(amplitudes[divisions - 1]);
}
#endif

// --- Utility routines -------------------------------------------------------

Colour WavelengthToXYZ(ClrReal lambda)
{
    Colour c;

    // calculate XYZ
    c[0] = kCIE_xbar.Sample(lambda);
    c[1] = kCIE_ybar.Sample(lambda);
    c[2] = kCIE_zbar.Sample(lambda);
    
    return(c);
}

// chromaticity stuff originally from utah's sunsky code

// 300-830 10nm
static Float tS0Amplitudes[54] =
{
    0.04,6.0,29.6,55.3,57.3,
    61.8,61.5,68.8,63.4,65.8,
    94.8,104.8,105.9,96.8,113.9,
    125.6,125.5,121.3,121.3,113.5,
    113.1,110.8,106.5,108.8,105.3,
    104.4,100.0,96.0,95.1,89.1,
    90.5,90.3,88.4,84.0,85.1,
    81.9,82.6,84.9,81.3,71.9,
    74.3,76.4,63.3,71.7,77.0,
    65.2,47.7,68.6,65.0,66.0,
    61.0,53.3,58.9,61.9
};

static Float tS1Amplitudes[54] =
{
    0.02,4.5,22.4,42.0,40.6,
    41.6,38.0,42.4,38.5,35.0,
    43.4,46.3,43.9,37.1,36.7,
    35.9,32.6,27.9,24.3,20.1,
    16.2,13.2,8.6,6.1,4.2,
    1.9,0.0,-1.6,-3.5,-3.5,
    -5.8,-7.2,-8.6,-9.5,-10.9,
    -10.7,-12.0,-14.0,-13.6,-12.0,
    -13.3,-12.9,-10.6,-11.6,-12.2,
    -10.2,-7.8,-11.2,-10.4,-10.6,
    -9.7,-8.3,-9.3,-9.8
};

static Float tS2Amplitudes[54] =
{
    0.0,2.0,4.0,8.5,7.8,
    6.7,5.3,6.1,3.0,1.2,
    -1.1,-0.5,-0.7,-1.2,-2.6,
    -2.9,-2.8,-2.6,-2.6,-1.8,
    -1.5,-1.3,-1.2,-1.0,-0.5,
    -0.3,0.0,0.2,0.5,2.1,
    3.2,4.1,4.7,5.1,6.7,
    7.3,8.6,9.8,10.2,8.3,
    9.6,8.5,7.0,7.6,8.0,
    6.7,5.2,7.4,6.8,7.0,
    6.4,5.5,6.1,6.5
};

RegularSpectralCurve    kS0Spectrum(tS0Amplitudes, 300, 830, 54);
RegularSpectralCurve    kS1Spectrum(tS1Amplitudes, 300, 830, 54);
RegularSpectralCurve    kS2Spectrum(tS2Amplitudes, 300, 830, 54);

ChromaticitySpectrum::ChromaticitySpectrum(ClrReal x, ClrReal y)
{
    M1 = (-1.3515 - 1.7703 * x +  5.9114 * y) / (0.0241 + 0.2562 * x - 0.7341 * y);
    M2 = ( 0.03   -31.4424 * x + 30.0717 * y) / (0.0241 + 0.2562 * x - 0.7341 * y);
}

Float ChromaticitySpectrum::Sample(Float lambda) const
{
    return(kS0Spectrum.Sample(lambda) + M1 * kS1Spectrum.Sample(lambda)
            + M2 * kS2Spectrum.Sample(lambda));
}

static Coord WavelengthToChroma(Float lambda)
{
    Colour c;

    // calculate XYZ
    c[0] = kCIE_xbar.Sample(lambda);
    c[1] = kCIE_ybar.Sample(lambda);
    c[2] = kCIE_zbar.Sample(lambda);

    c /= (c[0] + c[1] + c[2]);

    return(Coord(c[0], c[1]));
}

Void DrawChromaticityImage(Image &image)
// Plot CIE colour diagram
{
    Int     i, j;
    Colour  chromaClr;
    Coord   a1, a2, b1, b2;
    Float   aLambda, bLambda;
    Int     startX, endX;
    Bool    lastClip, clip;
    
    // the only tricky part to drawing the CIE
    // diagram is culling everything outside the (x,y)
    // curve defined by the pure wavelengths, and the
    // "line of purples" between the highest and lowest
    // wavelength in the spectrum.

    // to take care of this we use a1/a2/b1/b2 to step
    // through the spectral curve from either end, starting
    // from the beginning, and pick the correct start and
    // ending x values for each scanline.

    image.Clear(cGrey50);

    a1 = WavelengthToChroma(kCIE_Start);
    aLambda = kCIE_Start + kCIE_Delta;
    a2 = WavelengthToChroma(aLambda);
    b1 = a1;
    bLambda = kCIE_End;
    b2 = WavelengthToChroma(bLambda);
    
    for (j = 0; j < image.Height(); j++)
    { 
        Float y = (j + 0.5) / image.Height();
        
        if (y < a1[1])
            continue;
        // use 'while' here because some consecutive wavelength
        // samples will have the same (x,y) values
        while (y > a2[1])
        {
            a1 = a2;
            aLambda += kCIE_Delta;
            a2 = WavelengthToChroma(aLambda);
        }
        while (y > b2[1])
        {
            b1 = b2;
            bLambda -= kCIE_Delta;
            b2 = WavelengthToChroma(bLambda);
        }
    
        if (bLambda < aLambda)
            // all done
            break;

        startX = Int(image.Width() * mix(a1[0], a2[0], (y - a1[1]) / (a2[1] - a1[1])));
        endX = Int(image.Width() * mix(b1[0], b2[0], (y - b1[1]) / (b2[1] - b1[1])));
        if (endX >= image.Width())
            endX = image.Width() - 1;

        lastClip = true;
        
        for (i = startX; i < endX; i++)
        {
            Float x = (i + 0.5) / image.Width();

            chromaClr = Colour(x, y, (1.0 - x - y));
//          chromaClr = Colour(x / y, 1.0, (1.0 - x - y) / y); 
            clip = MinCmpt(csDisplay.fromXYZ * chromaClr) < 0.0;
            if (clip != lastClip)
            {
                image.SetPixel(i, j, cYellow);
                lastClip = clip;
            }
            else
            {
                chromaClr = csDisplay.XYZToGamutClip(chromaClr);
                image.SetPixel(i, j, chromaClr);
            }
        }
    }
}

Colour CIEXYZToLuv(const Colour &cXYZ, const Colour &cWhiteXYZ)
{
    ClrReal         lumRatio = cXYZ[1] / cWhiteXYZ[1];
    const Colour    denomRat(1.0, 15.0, 3.0);
    ClrReal         denom, u, v, lumStar;
    Colour          result;
    
    if (lumRatio <= 0.008856)
        lumStar = 903.3 * lumRatio;
    else
        lumStar = 116.0 * pow(lumRatio, 1.0 / 3.0) -  16.0;

    denom = dot(denomRat, cXYZ);
    u = 4.0 * cXYZ[0] / denom;
    v = 9.0 * cXYZ[1] / denom;
    denom = dot(denomRat, cWhiteXYZ);
    u -= 4.0 * cWhiteXYZ[0] / denom;
    v -= 9.0 * cWhiteXYZ[1] / denom;

    result[0] = lumStar;
    result[1] = 13.0 * lumStar * u;
    result[2] = 13.0 * lumStar * v;

    return(result);
}

BlackBodySpectrum::BlackBodySpectrum(Float temperature)
{
    temp = temperature;
}

Float BlackBodySpectrum::Sample(Float lambda) const
{
    Float   wavelength = lambda * 1e-9;
    Float   result;

    result = 3.74183e-16 * pow(wavelength, -5.0);
    result /= exp(1.4388e-2 / (wavelength * temp)) - 1.0;

    return(result);
}

---