5.2 The SampledSpectrum Class

SampledSpectrum uses the CoefficientSpectrum infrastructure to represent an SPD with uniformly spaced samples between a starting and an ending wavelength. The wavelength range covers from 400 nm to 700 nm—the range of the visual spectrum where the human visual system is most sensitive. The number of samples, 60, is generally more than enough to accurately represent complex SPDs for rendering. (See the “Further Reading” section for background on sampling rate requirements for SPDs.) Thus, the first sample represents the wavelength range left-bracket 400 comma 405 right-parenthesis , the second represents left-bracket 405 comma 410 right-parenthesis , and so forth. These values can easily be changed here as needed.

<<Spectrum Utility Declarations>>= 
static const int sampledLambdaStart = 400; static const int sampledLambdaEnd = 700; static const int nSpectralSamples = 60;

<<Spectrum Declarations>>+=  
class SampledSpectrum : public CoefficientSpectrum<nSpectralSamples> { public: <<SampledSpectrum Public Methods>> 
SampledSpectrum(Float v = 0.f) : CoefficientSpectrum(v) { } SampledSpectrum(const CoefficientSpectrum<nSpectralSamples> &v) : CoefficientSpectrum<nSpectralSamples>(v) { } static SampledSpectrum FromSampled(const Float *lambda, const Float *v, int n) { <<Sort samples if unordered, use sorted for returned spectrum>> 
if (!SpectrumSamplesSorted(lambda, v, n)) { std::vector<Float> slambda(&lambda[0], &lambda[n]); std::vector<Float> sv(&v[0], &v[n]); SortSpectrumSamples(&slambda[0], &sv[0], n); return FromSampled(&slambda[0], &sv[0], n); }
SampledSpectrum r; for (int i = 0; i < nSpectralSamples; ++i) { <<Compute average value of given SPD over i th sample’s range>> 
Float lambda0 = Lerp(Float(i) / Float(nSpectralSamples), sampledLambdaStart, sampledLambdaEnd); Float lambda1 = Lerp(Float(i + 1) / Float(nSpectralSamples), sampledLambdaStart, sampledLambdaEnd); r.c[i] = AverageSpectrumSamples(lambda, v, n, lambda0, lambda1);
} return r; } static void Init() { <<Compute XYZ matching functions for SampledSpectrum>> 
for (int i = 0; i < nSpectralSamples; ++i) { Float wl0 = Lerp(Float(i) / Float(nSpectralSamples), sampledLambdaStart, sampledLambdaEnd); Float wl1 = Lerp(Float(i + 1) / Float(nSpectralSamples), sampledLambdaStart, sampledLambdaEnd); X.c[i] = AverageSpectrumSamples(CIE_lambda, CIE_X, nCIESamples, wl0, wl1); Y.c[i] = AverageSpectrumSamples(CIE_lambda, CIE_Y, nCIESamples, wl0, wl1); Z.c[i] = AverageSpectrumSamples(CIE_lambda, CIE_Z, nCIESamples, wl0, wl1); }
<<Compute RGB to spectrum functions for SampledSpectrum>> 
for (int i = 0; i < nSpectralSamples; ++i) { Float wl0 = Lerp(Float(i) / Float(nSpectralSamples), sampledLambdaStart, sampledLambdaEnd); Float wl1 = Lerp(Float(i+1) / Float(nSpectralSamples), sampledLambdaStart, sampledLambdaEnd); rgbRefl2SpectWhite.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectWhite, nRGB2SpectSamples, wl0, wl1); rgbRefl2SpectCyan.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectCyan, nRGB2SpectSamples, wl0, wl1); rgbRefl2SpectMagenta.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectMagenta, nRGB2SpectSamples, wl0, wl1); rgbRefl2SpectYellow.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectYellow, nRGB2SpectSamples, wl0, wl1); rgbRefl2SpectRed.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectRed, nRGB2SpectSamples, wl0, wl1); rgbRefl2SpectGreen.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectGreen, nRGB2SpectSamples, wl0, wl1); rgbRefl2SpectBlue.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectBlue, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectWhite.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectWhite, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectCyan.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectCyan, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectMagenta.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectMagenta, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectYellow.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectYellow, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectRed.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectRed, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectGreen.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectGreen, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectBlue.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectBlue, nRGB2SpectSamples, wl0, wl1); }
} void ToXYZ(Float xyz[3]) const { xyz[0] = xyz[1] = xyz[2] = 0.f; for (int i = 0; i < nSpectralSamples; ++i) { xyz[0] += X.c[i] * c[i]; xyz[1] += Y.c[i] * c[i]; xyz[2] += Z.c[i] * c[i]; } Float scale = Float(sampledLambdaEnd - sampledLambdaStart) / Float(CIE_Y_integral * nSpectralSamples); xyz[0] *= scale; xyz[1] *= scale; xyz[2] *= scale; } Float y() const { Float yy = 0.f; for (int i = 0; i < nSpectralSamples; ++i) yy += Y.c[i] * c[i]; return yy * Float(sampledLambdaEnd - sampledLambdaStart) / Float(CIE_Y_integral * nSpectralSamples); } void ToRGB(Float rgb[3]) const { Float xyz[3]; ToXYZ(xyz); XYZToRGB(xyz, rgb); } RGBSpectrum ToRGBSpectrum() const; static SampledSpectrum FromRGB(const Float rgb[3], SpectrumType type = SpectrumType::Illuminant); static SampledSpectrum FromXYZ(const Float xyz[3], SpectrumType type = SpectrumType::Reflectance) { Float rgb[3]; XYZToRGB(xyz, rgb); return FromRGB(rgb, type); } SampledSpectrum(const RGBSpectrum &r, SpectrumType type = SpectrumType::Reflectance);
private: <<SampledSpectrum Private Data>> 
static SampledSpectrum X, Y, Z; static SampledSpectrum rgbRefl2SpectWhite, rgbRefl2SpectCyan; static SampledSpectrum rgbRefl2SpectMagenta, rgbRefl2SpectYellow; static SampledSpectrum rgbRefl2SpectRed, rgbRefl2SpectGreen; static SampledSpectrum rgbRefl2SpectBlue; static SampledSpectrum rgbIllum2SpectWhite, rgbIllum2SpectCyan; static SampledSpectrum rgbIllum2SpectMagenta, rgbIllum2SpectYellow; static SampledSpectrum rgbIllum2SpectRed, rgbIllum2SpectGreen; static SampledSpectrum rgbIllum2SpectBlue;
};

By inheriting from the CoefficientSpectrum class, SampledSpectrum automatically has all of the basic spectrum arithmetic operations defined earlier. The main methods left to define for it are those that initialize it from spectral data and convert the SPD it represents to other spectral representations (such as RGB). The constructor for initializing it with a constant SPD is straightforward.

<<SampledSpectrum Public Methods>>= 
SampledSpectrum(Float v = 0.f) : CoefficientSpectrum(v) { }

We will often be provided spectral data as a set of left-parenthesis lamda Subscript i Baseline comma v Subscript i Baseline right-parenthesis samples, where the i th sample has some value v Subscript i at wavelength lamda Subscript i . In general, the samples may have an irregular spacing and there may be more or fewer of them than the SampledSpectrum stores. (See the directory scenes/spds in the pbrt distribution for a variety of SPDs for use with pbrt, many of them with irregular sample spacing.)

The FromSampled() method takes arrays of SPD sample values v at given wavelengths lambda and uses them to define a piecewise linear function to represent the SPD. For each SPD sample in the SampledSpectrum, it uses the AverageSpectrumSamples() utility function, defined below, to compute the average of the piecewise linear function over the range of wavelengths that each SPD sample is responsible for.

<<SampledSpectrum Public Methods>>+=  
static SampledSpectrum FromSampled(const Float *lambda, const Float *v, int n) { <<Sort samples if unordered, use sorted for returned spectrum>> 
if (!SpectrumSamplesSorted(lambda, v, n)) { std::vector<Float> slambda(&lambda[0], &lambda[n]); std::vector<Float> sv(&v[0], &v[n]); SortSpectrumSamples(&slambda[0], &sv[0], n); return FromSampled(&slambda[0], &sv[0], n); }
SampledSpectrum r; for (int i = 0; i < nSpectralSamples; ++i) { <<Compute average value of given SPD over i th sample’s range>> 
Float lambda0 = Lerp(Float(i) / Float(nSpectralSamples), sampledLambdaStart, sampledLambdaEnd); Float lambda1 = Lerp(Float(i + 1) / Float(nSpectralSamples), sampledLambdaStart, sampledLambdaEnd); r.c[i] = AverageSpectrumSamples(lambda, v, n, lambda0, lambda1);
} return r; }

The AverageSpectrumSamples() function requires that the left-parenthesis lamda Subscript i Baseline comma v Subscript i Baseline right-parenthesis values be sorted by wavelength. The SpectrumSamplesSorted() function checks that they are; if not, SortSpectrumSamples() sorts them. Note that we allocate new storage for the sorted samples and do not modify the values passed in by the caller; in general, doing so would likely be unexpected behavior for a user of this function (who shouldn’t need to worry about these requirements of its specific implementation). We won’t include the implementations of either of these two functions here, as they are straightforward.

<<Sort samples if unordered, use sorted for returned spectrum>>= 
if (!SpectrumSamplesSorted(lambda, v, n)) { std::vector<Float> slambda(&lambda[0], &lambda[n]); std::vector<Float> sv(&v[0], &v[n]); SortSpectrumSamples(&slambda[0], &sv[0], n); return FromSampled(&slambda[0], &sv[0], n); }

In order to compute the value for the i th spectral sample, we compute the range of wavelengths that it’s responsible for—lambda0 to lambda1—and use the AverageSpectrumSamples() function to compute the average value of the given piecewise linear SPD over that range. This is a 1D instance of sampling and reconstruction, a topic that will be discussed in more detail in Chapter 7.

<<Compute average value of given SPD over i th sample’s range>>= 
Float lambda0 = Lerp(Float(i) / Float(nSpectralSamples), sampledLambdaStart, sampledLambdaEnd); Float lambda1 = Lerp(Float(i + 1) / Float(nSpectralSamples), sampledLambdaStart, sampledLambdaEnd); r.c[i] = AverageSpectrumSamples(lambda, v, n, lambda0, lambda1);

Figure 5.2 shows the basic approach taken by AverageSpectrumSamples(): it iterates over each of the linear segments between samples that are partially or fully within the range of wavelengths, lambdaStart to lambdaEnd. For each such segment, it computes the average value over its range, scales the average by the wavelength range the segment covers, and accumulates a sum of these values. The final average value is this sum divided by the total wavelength range.

Figure 5.2: When resampling an irregularly defined SPD, we need to compute the average value of the piecewise linear function defined by the SPD samples. Here, we want to average the value from 500 nm to 600 nm—the shaded region under the plot. The AverageSpectrumSamples() function does this by computing the areas of each of the regions denoted by dashed lines in this figure.

<<Spectrum Method Definitions>>= 
Float AverageSpectrumSamples(const Float *lambda, const Float *vals, int n, Float lambdaStart, Float lambdaEnd) { <<Handle cases with out-of-bounds range or single sample only>> 
if (lambdaEnd <= lambda[0]) return vals[0]; if (lambdaStart >= lambda[n - 1]) return vals[n - 1]; if (n == 1) return vals[0];
Float sum = 0; <<Add contributions of constant segments before/after samples>> 
if (lambdaStart < lambda[0]) sum += vals[0] * (lambda[0] - lambdaStart); if (lambdaEnd > lambda[n-1]) sum += vals[n - 1] * (lambdaEnd - lambda[n - 1]);
<<Advance to first relevant wavelength segment>> 
int i = 0; while (lambdaStart > lambda[i + 1]) ++i;
<<Loop over wavelength sample segments and add contributions>> 
auto interp = [lambda, vals](Float w, int i) { return Lerp((w - lambda[i]) / (lambda[i + 1] - lambda[i]), vals[i], vals[i + 1]); }; for (; i+1 < n && lambdaEnd >= lambda[i]; ++i) { Float segLambdaStart = std::max(lambdaStart, lambda[i]); Float segLambdaEnd = std::min(lambdaEnd, lambda[i + 1]); sum += 0.5 * (interp(segLambdaStart, i) + interp(segLambdaEnd, i)) * (segLambdaEnd - segLambdaStart); }
return sum / (lambdaEnd - lambdaStart); }

The function starts by checking for and handling the edge cases where the range of wavelengths to average over is outside the range of provided wavelengths or the case where there is only a single sample, in which case the average value is trivially computed. We assume that the SPD has a constant value (the values at the two respective endpoints) outside of the provided sample range; if this isn’t a reasonable assumption for a particular set of data, the data provided should have explicit values of (for example) 0 at the endpoints.

<<Handle cases with out-of-bounds range or single sample only>>= 
if (lambdaEnd <= lambda[0]) return vals[0]; if (lambdaStart >= lambda[n - 1]) return vals[n - 1]; if (n == 1) return vals[0];

Having handled these cases, the next step is to check to see if part of the range to average over goes beyond the first and/or last sample value. If so, we accumulate the contribution of the constant segment(s), scaled by the out-of-bounds wavelength range.

<<Add contributions of constant segments before/after samples>>= 
if (lambdaStart < lambda[0]) sum += vals[0] * (lambda[0] - lambdaStart); if (lambdaEnd > lambda[n-1]) sum += vals[n - 1] * (lambdaEnd - lambda[n - 1]);

And now we advance to the first index i where the starting wavelength of the interpolation range overlaps the segment from lamda Subscript i to lamda Subscript i plus 1 . A more efficient implementation would use a binary search rather than a linear search here. However, this code is currently only called at scene initialization time, so the lack of these optimizations doesn’t currently impact rendering performance.

<<Advance to first relevant wavelength segment>>= 
int i = 0; while (lambdaStart > lambda[i + 1]) ++i;

The loop below iterates over each of the linear segments that the averaging range overlaps. For each one, it computes the average value over the wavelength range segLambdaStart to segLambdaEnd by averaging the values of the function at those two points. The values in turn are computed by interp(), a lambda function that linearly interpolates between the two endpoints at the given wavelength.

The std::min() and std::max() calls below compute the wavelength range to average over within the segment; note that they naturally handle the cases where lambdaStart, lambdaEnd, or both of them are within the current segment.

<<Loop over wavelength sample segments and add contributions>>= 
auto interp = [lambda, vals](Float w, int i) { return Lerp((w - lambda[i]) / (lambda[i + 1] - lambda[i]), vals[i], vals[i + 1]); }; for (; i+1 < n && lambdaEnd >= lambda[i]; ++i) { Float segLambdaStart = std::max(lambdaStart, lambda[i]); Float segLambdaEnd = std::min(lambdaEnd, lambda[i + 1]); sum += 0.5 * (interp(segLambdaStart, i) + interp(segLambdaEnd, i)) * (segLambdaEnd - segLambdaStart); }

5.2.1 XYZ Color

A remarkable property of the human visual system makes it possible to represent colors for human perception with just three floating-point numbers. The tristimulus theory of color perception says that all visible SPDs can be accurately represented for human observers with three values, x Subscript lamda , y Subscript lamda , and z Subscript lamda . Given an SPD upper S left-parenthesis lamda right-parenthesis , these values are computed by integrating its product with the spectral matching curves upper X left-parenthesis lamda right-parenthesis , upper Y left-parenthesis lamda right-parenthesis , and upper Z left-parenthesis lamda right-parenthesis :

StartLayout 1st Row 1st Column x Subscript lamda 2nd Column equals StartFraction 1 Over integral upper Y left-parenthesis lamda right-parenthesis normal d lamda Subscript Baseline EndFraction integral Underscript lamda Endscripts upper S left-parenthesis lamda right-parenthesis upper X left-parenthesis lamda right-parenthesis normal d lamda Subscript Baseline 2nd Row 1st Column y Subscript lamda 2nd Column equals StartFraction 1 Over integral upper Y left-parenthesis lamda right-parenthesis normal d lamda Subscript Baseline EndFraction integral Underscript lamda Endscripts upper S left-parenthesis lamda right-parenthesis upper Y left-parenthesis lamda right-parenthesis normal d lamda Subscript Baseline 3rd Row 1st Column z Subscript lamda 2nd Column equals StartFraction 1 Over integral upper Y left-parenthesis lamda right-parenthesis normal d lamda Subscript Baseline EndFraction integral Underscript lamda Endscripts upper S left-parenthesis lamda right-parenthesis upper Z left-parenthesis lamda right-parenthesis normal d lamda Subscript Baseline period EndLayout
(5.1)

These curves were determined by the Commission Internationale de l’Éclairage (CIE) standards body after a series of experiments with human test subjects and are graphed in Figure 5.3. It is believed that these matching curves are generally similar to the responses of the three types of color-sensitive cones in the human retina. Remarkably, SPDs with substantially different distributions may have very similar x Subscript lamda , y Subscript lamda , and z Subscript lamda values. To the human observer, such SPDs actually appear the same visually. Pairs of such spectra are called metamers.

Figure 5.3: Computing the XYZ Values for an Arbitrary SPD. The SPD is multiplied by each of the three matching curves and integrated over their non-zero extent to compute the values x Subscript lamda , y Subscript lamda , and z Subscript lamda , using Equation (5.1).

This brings us to a subtle point about representations of spectral power distributions. Most color spaces attempt to model colors that are visible to humans and therefore use only three coefficients, exploiting the tristimulus theory of color perception. Although XYZ works well to represent a given SPD to be displayed for a human observer, it is not a particularly good set of basis functions for spectral computation. For example, although XYZ values would work well to describe the perceived color of lemon skin or a fluorescent light individually (recall Figure 5.1), the product of their respective XYZ values is likely to give a noticeably different XYZ color than the XYZ value computed by multiplying more accurate representations of their SPDs and then computing the XYZ value.

pbrt provides the values of the standard upper X left-parenthesis lamda right-parenthesis , upper Y left-parenthesis lamda right-parenthesis , and upper Z left-parenthesis lamda right-parenthesis response curves sampled at 1-nm increments from 360 nm to 830 nm. The wavelengths of the i th sample in the arrays below are given by the i th element of CIE_lambda; having the wavelengths of the samples explicitly represented in this way makes it easy to pass the XYZ samples into functions like AverageSpectrumSamples() that take such an array of wavelengths as a parameter.

<<Spectral Data Declarations>>= 
static const int nCIESamples = 471; extern const Float CIE_X[nCIESamples]; extern const Float CIE_Y[nCIESamples]; extern const Float CIE_Z[nCIESamples]; extern const Float CIE_lambda[nCIESamples];

SampledSpectrum uses these samples to compute the XYZ matching curves in its spectral representation (i.e., themselves as SampledSpectrums).

<<SampledSpectrum Private Data>>= 
static SampledSpectrum X, Y, Z;

The SampledSpectrum XYZ matching curves are computed in the SampledSpectrum::Init() method, which is called at system startup time by the pbrtInit() function defined in Section A.2.

<<SampledSpectrum Public Methods>>+=  
static void Init() { <<Compute XYZ matching functions for SampledSpectrum>> 
for (int i = 0; i < nSpectralSamples; ++i) { Float wl0 = Lerp(Float(i) / Float(nSpectralSamples), sampledLambdaStart, sampledLambdaEnd); Float wl1 = Lerp(Float(i + 1) / Float(nSpectralSamples), sampledLambdaStart, sampledLambdaEnd); X.c[i] = AverageSpectrumSamples(CIE_lambda, CIE_X, nCIESamples, wl0, wl1); Y.c[i] = AverageSpectrumSamples(CIE_lambda, CIE_Y, nCIESamples, wl0, wl1); Z.c[i] = AverageSpectrumSamples(CIE_lambda, CIE_Z, nCIESamples, wl0, wl1); }
<<Compute RGB to spectrum functions for SampledSpectrum>> 
for (int i = 0; i < nSpectralSamples; ++i) { Float wl0 = Lerp(Float(i) / Float(nSpectralSamples), sampledLambdaStart, sampledLambdaEnd); Float wl1 = Lerp(Float(i+1) / Float(nSpectralSamples), sampledLambdaStart, sampledLambdaEnd); rgbRefl2SpectWhite.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectWhite, nRGB2SpectSamples, wl0, wl1); rgbRefl2SpectCyan.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectCyan, nRGB2SpectSamples, wl0, wl1); rgbRefl2SpectMagenta.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectMagenta, nRGB2SpectSamples, wl0, wl1); rgbRefl2SpectYellow.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectYellow, nRGB2SpectSamples, wl0, wl1); rgbRefl2SpectRed.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectRed, nRGB2SpectSamples, wl0, wl1); rgbRefl2SpectGreen.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectGreen, nRGB2SpectSamples, wl0, wl1); rgbRefl2SpectBlue.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBRefl2SpectBlue, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectWhite.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectWhite, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectCyan.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectCyan, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectMagenta.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectMagenta, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectYellow.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectYellow, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectRed.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectRed, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectGreen.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectGreen, nRGB2SpectSamples, wl0, wl1); rgbIllum2SpectBlue.c[i] = AverageSpectrumSamples(RGB2SpectLambda, RGBIllum2SpectBlue, nRGB2SpectSamples, wl0, wl1); }
}

<<General pbrt Initialization>>= 

Given the wavelength range and number of samples for SampledSpectrum, computing the values of the matching functions for each sample is just a matter of computing the sample’s wavelength range and using the AverageSpectrumSamples() routine.

<<Compute XYZ matching functions for SampledSpectrum>>= 
for (int i = 0; i < nSpectralSamples; ++i) { Float wl0 = Lerp(Float(i) / Float(nSpectralSamples), sampledLambdaStart, sampledLambdaEnd); Float wl1 = Lerp(Float(i + 1) / Float(nSpectralSamples), sampledLambdaStart, sampledLambdaEnd); X.c[i] = AverageSpectrumSamples(CIE_lambda, CIE_X, nCIESamples, wl0, wl1); Y.c[i] = AverageSpectrumSamples(CIE_lambda, CIE_Y, nCIESamples, wl0, wl1); Z.c[i] = AverageSpectrumSamples(CIE_lambda, CIE_Z, nCIESamples, wl0, wl1); }

All Spectrum implementations in pbrt must provide a method that converts their SPD to left-parenthesis x Subscript lamda Baseline comma y Subscript lamda Baseline comma z Subscript lamda Baseline right-parenthesis coefficients. This method is called, for example, in the process of updating pixels in the image. When a Spectrum representing the light carried along a ray from the camera is provided to the Film, the Film converts the SPD into XYZ coefficients as a first step in the process of finally turning them into RGB values used for storage and/or display.

To compute XYZ coefficients, SampledSpectrum computes the integrals from Equation (5.1) with a Riemann sum:

x Subscript lamda Baseline almost-equals StartFraction 1 Over integral upper Y left-parenthesis lamda right-parenthesis normal d lamda Subscript Baseline EndFraction StartFraction lamda Subscript normal e normal n normal d Baseline minus lamda Subscript normal s normal t normal a normal r normal t Baseline Over upper N EndFraction sigma-summation Underscript i equals 0 Overscript upper N minus 1 Endscripts upper X Subscript i Baseline c Subscript i Baseline comma

and so forth.

<<SampledSpectrum Public Methods>>+=  
void ToXYZ(Float xyz[3]) const { xyz[0] = xyz[1] = xyz[2] = 0.f; for (int i = 0; i < nSpectralSamples; ++i) { xyz[0] += X.c[i] * c[i]; xyz[1] += Y.c[i] * c[i]; xyz[2] += Z.c[i] * c[i]; } Float scale = Float(sampledLambdaEnd - sampledLambdaStart) / Float(CIE_Y_integral * nSpectralSamples); xyz[0] *= scale; xyz[1] *= scale; xyz[2] *= scale; }

The value of the integral integral upper Y left-parenthesis lamda right-parenthesis normal d lamda Subscript is used in a number of calculations like these; it’s therefore useful to have its value available directly via the CIE_Y_integral constant.

<<Spectral Data Declarations>>+=  
static const Float CIE_Y_integral = 106.856895;

The y coefficient of XYZ color is closely related to luminance, which measures the perceived brightness of a color. Luminance is described in more detail in Section 5.4.3. We provide a method to compute y alone in a separate method as often only the luminance of a spectrum is desired. (For example, some of the light transport algorithms in Chapters 14, 15, and 16 use luminance as a measure of the relative importance of light-carrying paths through the scene.)

<<SampledSpectrum Public Methods>>+=  
Float y() const { Float yy = 0.f; for (int i = 0; i < nSpectralSamples; ++i) yy += Y.c[i] * c[i]; return yy * Float(sampledLambdaEnd - sampledLambdaStart) / Float(CIE_Y_integral * nSpectralSamples); }

5.2.2 RGB Color

When we display an RGB color on a display, the spectrum that is actually displayed is basically determined by the weighted sum of three spectral response curves, one for each of red, green, and blue, as emitted by the display’s phosphors, LED or LCD elements, or plasma cells. Figure 5.4 plots the red, green, and blue distributions emitted by a LED display and a LCD display; note that they are remarkably different.

Figure 5.4: Red, Green, and Blue Emission Curves for an LCD Display and an LED Display. The first plot shows the curves for an LCD display, and the second shows them for an LED. These two displays have quite different emission profiles. (Data courtesy of X-Rite, Inc.)

Figure 5.5 in turn shows the SPDs that result from displaying the RGB color left-parenthesis 0.6 comma 0.3 comma 0.2 right-parenthesis on those displays. Not surprisingly, the resulting SPDs are quite different as well. This example illustrates that using RGB values provided by the user to describe a particular color is actually only meaningful given knowledge of the characteristics of the display they were using when they selected the RGB values.

Figure 5.5: SPDs from Displaying the RGB Color left-parenthesis 0.6 comma 0.3 comma 0.2 right-parenthesis on LED and LCD Displays. The resulting emitted SPDs are remarkably different, even given the same RGB values, due to the different emission curves illustrated in Figure 5.4.

Given an left-parenthesis x Subscript lamda Baseline comma y Subscript lamda Baseline comma z Subscript lamda Baseline right-parenthesis representation of an SPD, we can convert it to corresponding RGB coefficients, given the choice of a particular set of SPDs that define red, green, and blue for a display of interest. Given the spectral response curves upper R left-parenthesis lamda right-parenthesis , upper G left-parenthesis lamda right-parenthesis , and upper B left-parenthesis lamda right-parenthesis , for a particular display, RGB coefficients can be computed by integrating the response curves with the SPD upper S left-parenthesis lamda right-parenthesis and using the tristimulus theory of color perception:

StartLayout 1st Row 1st Column r 2nd Column equals integral upper R left-parenthesis lamda right-parenthesis upper S left-parenthesis lamda right-parenthesis normal d lamda Subscript Baseline equals integral upper R left-parenthesis lamda right-parenthesis left-parenthesis x Subscript lamda Baseline upper X left-parenthesis lamda right-parenthesis plus y Subscript lamda Baseline upper Y left-parenthesis lamda right-parenthesis plus z Subscript lamda Baseline upper Z left-parenthesis lamda right-parenthesis right-parenthesis normal d lamda Subscript Baseline 2nd Row 1st Column Blank 2nd Column equals x Subscript lamda Baseline integral upper R left-parenthesis lamda right-parenthesis upper X left-parenthesis lamda right-parenthesis normal d lamda Subscript Baseline plus y Subscript lamda Baseline integral upper R left-parenthesis lamda right-parenthesis upper Y left-parenthesis lamda right-parenthesis normal d lamda Subscript Baseline plus z Subscript lamda Baseline integral upper R left-parenthesis lamda right-parenthesis upper Z left-parenthesis lamda right-parenthesis normal d lamda Subscript Baseline period EndLayout

The integrals of the products of upper R left-parenthesis lamda right-parenthesis upper X left-parenthesis lamda right-parenthesis and so forth can be precomputed for given response curves, making it possible to express the full conversion as a matrix:

Start 3 By 1 Matrix 1st Row r 2nd Row g 3rd Row b EndMatrix equals Start 3 By 3 Matrix 1st Row 1st Column integral upper R left-parenthesis lamda right-parenthesis upper X left-parenthesis lamda right-parenthesis normal d lamda Subscript Baseline 2nd Column integral upper R left-parenthesis lamda right-parenthesis upper Y left-parenthesis lamda right-parenthesis normal d lamda Subscript Baseline 3rd Column integral upper R left-parenthesis lamda right-parenthesis upper Z left-parenthesis lamda right-parenthesis normal d lamda Subscript Baseline 2nd Row 1st Column integral upper G left-parenthesis lamda right-parenthesis upper X left-parenthesis lamda right-parenthesis normal d lamda Subscript Baseline 2nd Column integral upper G left-parenthesis lamda right-parenthesis upper Y left-parenthesis lamda right-parenthesis normal d lamda Subscript Baseline 3rd Column integral upper G left-parenthesis lamda right-parenthesis upper Z left-parenthesis lamda right-parenthesis normal d lamda Subscript Baseline 3rd Row 1st Column integral upper B left-parenthesis lamda right-parenthesis upper X left-parenthesis lamda right-parenthesis normal d lamda Subscript Baseline 2nd Column integral upper B left-parenthesis lamda right-parenthesis upper Y left-parenthesis lamda right-parenthesis normal d lamda Subscript Baseline 3rd Column integral upper B left-parenthesis lamda right-parenthesis upper Z left-parenthesis lamda right-parenthesis normal d lamda Subscript Baseline EndMatrix Start 3 By 1 Matrix 1st Row x Subscript lamda Baseline 2nd Row y Subscript lamda Baseline 3rd Row z Subscript lamda Baseline EndMatrix period

The conversion routines implemented in pbrt are based on a standard set of these RGB spectra that has been defined for high-definition television.

<<Spectrum Utility Declarations>>+=  
inline void XYZToRGB(const Float xyz[3], Float rgb[3]) { rgb[0] = 3.240479f*xyz[0] - 1.537150f*xyz[1] - 0.498535f*xyz[2]; rgb[1] = -0.969256f*xyz[0] + 1.875991f*xyz[1] + 0.041556f*xyz[2]; rgb[2] = 0.055648f*xyz[0] - 0.204043f*xyz[1] + 1.057311f*xyz[2]; }

The inverse of this matrix gives coefficients to convert given RGB values expressed with respect to a particular set of RGB response curves to left-parenthesis x Subscript lamda Baseline comma y Subscript lamda Baseline comma z Subscript lamda Baseline right-parenthesis coefficients.

<<Spectrum Utility Declarations>>+=  
inline void RGBToXYZ(const Float rgb[3], Float xyz[3]) { xyz[0] = 0.412453f*rgb[0] + 0.357580f*rgb[1] + 0.180423f*rgb[2]; xyz[1] = 0.212671f*rgb[0] + 0.715160f*rgb[1] + 0.072169f*rgb[2]; xyz[2] = 0.019334f*rgb[0] + 0.119193f*rgb[1] + 0.950227f*rgb[2]; }

Given these functions, a SampledSpectrum can convert to RGB coefficients by first converting to XYZ and then using the XYZToRGB() utility function.

<<SampledSpectrum Public Methods>>+=  
void ToRGB(Float rgb[3]) const { Float xyz[3]; ToXYZ(xyz); XYZToRGB(xyz, rgb); }

An RGBSpectrum can also be created easily, using the ToRGBSpectrum() method.

<<SampledSpectrum Public Methods>>+=  
RGBSpectrum ToRGBSpectrum() const;

Going the other way and converting from RGB or XYZ values to a SPD isn’t as easy: the problem is highly under-constrained. Recall that an infinite number of different SPDs have the same left-parenthesis x Subscript lamda Baseline comma y Subscript lamda Baseline comma z Subscript lamda Baseline right-parenthesis (and thus RGB) coefficients. Thus, given an RGB or left-parenthesis x Subscript lamda Baseline comma y Subscript lamda Baseline comma z Subscript lamda Baseline right-parenthesis value, there are an infinite number of possible SPDs that could be chosen for it. There are a number of desirable criteria that we’d like a conversion function to have:

  • If all of the RGB coefficients have the same value, the resulting SPD should be constant.
  • In general, it’s desirable that the computed SPD be smooth. Most real-world objects have relatively smooth spectra. (The main source of spiky spectra is light sources, especially fluorescents. Fortunately, actual spectral data are more commonly available for illuminants than they are for reflectances.)

The smoothness goal is one of the reasons why constructing an SPD as a weighted sum of a display’s upper R left-parenthesis lamda right-parenthesis , upper G left-parenthesis lamda right-parenthesis , and upper B left-parenthesis lamda right-parenthesis SPDs is not a good solution: as shown in Figure 5.4, those functions are generally irregular and spiky, and a weighted sum of them will thus not be a very smooth SPD. Although the result will be a metamer of the given RGB values, it’s likely not an accurate representation of the SPD of the actual object.

Here we implement a method for converting RGBs to SPDs suggested by Smits (1999) that tries to achieve the goals above. This approach is based on the observation that a good start is to compute individual SPDs for red, green, and blue that are smooth and such that computing the weighted sum of them with the given RGB coefficients and then converting back to RGB give a result that is close to the original RGB coefficients. He found such spectra through a numerical optimization procedure.

Smits observed that two additional improvements could be made to this basic approach. First, rather than representing constant spectra by the sums of the computed red, green, and blue SPDs, the sum of which isn’t exactly constant, it’s better to represent constant spectra with constant SPDs. Second, mixtures of colors like yellow (a mixture of red and green) that are a mixture of two of the primaries are better represented by their own precomputed smooth SPDs rather than the sum of SPDs for the two corresponding primaries.

The following arrays store SPDs that meet these criteria, with their samples’ wavelengths in RGB2SpectLambda[] (these data were generated by Karl vom Berge).

<<Spectral Data Declarations>>+=  
static const int nRGB2SpectSamples = 32; extern const Float RGB2SpectLambda[nRGB2SpectSamples]; extern const Float RGBRefl2SpectWhite[nRGB2SpectSamples]; extern const Float RGBRefl2SpectCyan[nRGB2SpectSamples]; extern const Float RGBRefl2SpectMagenta[nRGB2SpectSamples]; extern const Float RGBRefl2SpectYellow[nRGB2SpectSamples]; extern const Float RGBRefl2SpectRed[nRGB2SpectSamples]; extern const Float RGBRefl2SpectGreen[nRGB2SpectSamples]; extern const Float RGBRefl2SpectBlue[nRGB2SpectSamples];

If a given RGB color describes illumination from a light source, better results are achieved if the conversion tables are computed using the spectral power distribution of a representative illumination source to define “white” rather than using a constant spectrum as they are for the tables above that are used for reflectances. The RGBIllum2Spect arrays use the D65 spectral power distribution, which has been standardized by the CIE to represent midday sunlight. (The D65 illuminant will be discussed more in Section 12.1.2.)

<<Spectral Data Declarations>>+= 
extern const Float RGBIllum2SpectWhite[nRGB2SpectSamples]; extern const Float RGBIllum2SpectCyan[nRGB2SpectSamples]; extern const Float RGBIllum2SpectMagenta[nRGB2SpectSamples]; extern const Float RGBIllum2SpectYellow[nRGB2SpectSamples]; extern const Float RGBIllum2SpectRed[nRGB2SpectSamples]; extern const Float RGBIllum2SpectGreen[nRGB2SpectSamples]; extern const Float RGBIllum2SpectBlue[nRGB2SpectSamples];

The fragment <<Compute RGB to spectrum functions for SampledSpectrum>>, which is called from SampledSpectrum::Init(), isn’t included here; it initializes the following SampledSpectrum values by resampling the RGBRefl2Spect and RGBIllum2Spect distributions using the AverageSpectrumSamples() function.

<<SampledSpectrum Private Data>>+=  
static SampledSpectrum rgbRefl2SpectWhite, rgbRefl2SpectCyan; static SampledSpectrum rgbRefl2SpectMagenta, rgbRefl2SpectYellow; static SampledSpectrum rgbRefl2SpectRed, rgbRefl2SpectGreen; static SampledSpectrum rgbRefl2SpectBlue;

<<SampledSpectrum Private Data>>+= 
static SampledSpectrum rgbIllum2SpectWhite, rgbIllum2SpectCyan; static SampledSpectrum rgbIllum2SpectMagenta, rgbIllum2SpectYellow; static SampledSpectrum rgbIllum2SpectRed, rgbIllum2SpectGreen; static SampledSpectrum rgbIllum2SpectBlue;

The SampledSpectrum::FromRGB() method converts from the given RGB values to a full SPD. In addition to the RGB values, it takes an enumeration value that denotes whether the RGB value represents surface reflectance or an illuminant; the corresponding rgbIllum2Spect or rgbRefl2Spect values are used for the conversion.

<<Spectrum Utility Declarations>>+= 
enum class SpectrumType { Reflectance, Illuminant };

<<Spectrum Method Definitions>>+=  
SampledSpectrum SampledSpectrum::FromRGB(const Float rgb[3], SpectrumType type) { SampledSpectrum r; if (type == SpectrumType::Reflectance) { <<Convert reflectance spectrum to RGB>> 
if (rgb[0] <= rgb[1] && rgb[0] <= rgb[2]) { <<Compute reflectance SampledSpectrum with rgb[0] as minimum>> 
r += rgb[0] * rgbRefl2SpectWhite; if (rgb[1] <= rgb[2]) { r += (rgb[1] - rgb[0]) * rgbRefl2SpectCyan; r += (rgb[2] - rgb[1]) * rgbRefl2SpectBlue; } else { r += (rgb[2] - rgb[0]) * rgbRefl2SpectCyan; r += (rgb[1] - rgb[2]) * rgbRefl2SpectGreen; }
} else if (rgb[1] <= rgb[0] && rgb[1] <= rgb[2]) { <<Compute reflectance SampledSpectrum with rgb[1] as minimum>> 
r += rgb[1] * rgbRefl2SpectWhite; if (rgb[0] <= rgb[2]) { r += (rgb[0] - rgb[1]) * rgbRefl2SpectMagenta; r += (rgb[2] - rgb[0]) * rgbRefl2SpectBlue; } else { r += (rgb[2] - rgb[1]) * rgbRefl2SpectMagenta; r += (rgb[0] - rgb[2]) * rgbRefl2SpectRed; }
} else { <<Compute reflectance SampledSpectrum with rgb[2] as minimum>> 
r += rgb[2] * rgbRefl2SpectWhite; if (rgb[0] <= rgb[1]) { r += (rgb[0] - rgb[2]) * rgbRefl2SpectYellow; r += (rgb[1] - rgb[0]) * rgbRefl2SpectGreen; } else { r += (rgb[1] - rgb[2]) * rgbRefl2SpectYellow; r += (rgb[0] - rgb[1]) * rgbRefl2SpectRed; }
}
} else { <<Convert illuminant spectrum to RGB>> 
if (rgb[0] <= rgb[1] && rgb[0] <= rgb[2]) { <<Compute illuminant SampledSpectrum with rgb[0] as minimum>> 
r += rgb[0] * rgbIllum2SpectWhite; if (rgb[1] <= rgb[2]) { r += (rgb[1] - rgb[0]) * rgbIllum2SpectCyan; r += (rgb[2] - rgb[1]) * rgbIllum2SpectBlue; } else { r += (rgb[2] - rgb[0]) * rgbIllum2SpectCyan; r += (rgb[1] - rgb[2]) * rgbIllum2SpectGreen; }
} else if (rgb[1] <= rgb[0] && rgb[1] <= rgb[2]) { <<Compute illuminant SampledSpectrum with rgb[1] as minimum>> 
r += rgb[1] * rgbIllum2SpectWhite; if (rgb[0] <= rgb[2]) { r += (rgb[0] - rgb[1]) * rgbIllum2SpectMagenta; r += (rgb[2] - rgb[0]) * rgbIllum2SpectBlue; } else { r += (rgb[2] - rgb[1]) * rgbIllum2SpectMagenta; r += (rgb[0] - rgb[2]) * rgbIllum2SpectRed; }
} else { <<Compute illuminant SampledSpectrum with rgb[2] as minimum>> 
r += rgb[2] * rgbIllum2SpectWhite; if (rgb[0] <= rgb[1]) { r += (rgb[0] - rgb[2]) * rgbIllum2SpectYellow; r += (rgb[1] - rgb[0]) * rgbIllum2SpectGreen; } else { r += (rgb[1] - rgb[2]) * rgbIllum2SpectYellow; r += (rgb[0] - rgb[1]) * rgbIllum2SpectRed; }
}
} return r.Clamp(); }

Here we’ll show the conversion process for reflectances. The computation for illuminants is the same, just using the different conversion values. First, the implementation determines whether the red, green, or blue channel is the smallest.

<<Convert reflectance spectrum to RGB>>= 
if (rgb[0] <= rgb[1] && rgb[0] <= rgb[2]) { <<Compute reflectance SampledSpectrum with rgb[0] as minimum>> 
r += rgb[0] * rgbRefl2SpectWhite; if (rgb[1] <= rgb[2]) { r += (rgb[1] - rgb[0]) * rgbRefl2SpectCyan; r += (rgb[2] - rgb[1]) * rgbRefl2SpectBlue; } else { r += (rgb[2] - rgb[0]) * rgbRefl2SpectCyan; r += (rgb[1] - rgb[2]) * rgbRefl2SpectGreen; }
} else if (rgb[1] <= rgb[0] && rgb[1] <= rgb[2]) { <<Compute reflectance SampledSpectrum with rgb[1] as minimum>> 
r += rgb[1] * rgbRefl2SpectWhite; if (rgb[0] <= rgb[2]) { r += (rgb[0] - rgb[1]) * rgbRefl2SpectMagenta; r += (rgb[2] - rgb[0]) * rgbRefl2SpectBlue; } else { r += (rgb[2] - rgb[1]) * rgbRefl2SpectMagenta; r += (rgb[0] - rgb[2]) * rgbRefl2SpectRed; }
} else { <<Compute reflectance SampledSpectrum with rgb[2] as minimum>> 
r += rgb[2] * rgbRefl2SpectWhite; if (rgb[0] <= rgb[1]) { r += (rgb[0] - rgb[2]) * rgbRefl2SpectYellow; r += (rgb[1] - rgb[0]) * rgbRefl2SpectGreen; } else { r += (rgb[1] - rgb[2]) * rgbRefl2SpectYellow; r += (rgb[0] - rgb[1]) * rgbRefl2SpectRed; }
}

Here is the code for the case of a red component being the smallest. (The cases for green and blue are analogous and not included in the book here.) If red is the smallest, we know that green and blue have greater values than red. As such, we can start to convert the final SPD to return by assigning to it the value of the red component times the white spectrum in rgbRefl2SpectWhite. Having done this, the remaining RGB value left to process is left-parenthesis 0 comma g minus r comma b minus r right-parenthesis . The code in turn determines which of the remaining two components is the smallest. This value, times the cyan (green and blue) spectrum, is added to the result and we’re left with either left-parenthesis 0 comma g minus b comma 0 right-parenthesis or left-parenthesis 0 comma 0 comma b minus g right-parenthesis . Based on whether the green or blue channel is non-zero, the green or blue SPD is scaled by the remainder and the conversion is complete.

<<Compute reflectance SampledSpectrum with rgb[0] as minimum>>= 
r += rgb[0] * rgbRefl2SpectWhite; if (rgb[1] <= rgb[2]) { r += (rgb[1] - rgb[0]) * rgbRefl2SpectCyan; r += (rgb[2] - rgb[1]) * rgbRefl2SpectBlue; } else { r += (rgb[2] - rgb[0]) * rgbRefl2SpectCyan; r += (rgb[1] - rgb[2]) * rgbRefl2SpectGreen; }

Given the method to convert from RGB, converting from XYZ color is easy. We first convert from XYZ to RGB and then use the FromRGB() method.

<<SampledSpectrum Public Methods>>+= 
static SampledSpectrum FromXYZ(const Float xyz[3], SpectrumType type = SpectrumType::Reflectance) { Float rgb[3]; XYZToRGB(xyz, rgb); return FromRGB(rgb, type); }

Finally, we provide a constructor that converts from an instance of the RGBSpectrum class, again using the infrastructure above.

<<Spectrum Method Definitions>>+=  
SampledSpectrum::SampledSpectrum(const RGBSpectrum &r, SpectrumType t) { Float rgb[3]; r.ToRGB(rgb); *this = SampledSpectrum::FromRGB(rgb, t); }