## 4.5 Representing Spectral Distributions

Spectral distributions in the real world can be complex;
we have already seen a variety of complex emission spectra and
Figure 4.16 shows a graph of the spectral distribution of
the reflectance of lemon skin. In order to render images of scenes that include a variety
of complex spectra, a renderer must have efficient and accurate
representations of spectral distributions. This section will introduce
`pbrt`’s abstractions for representing and performing computation with them;
the corresponding code can be found in the files `util/spectrum.h` and
`util/spectrum.cpp`.

We will start by defining constants that give the range of visible
wavelengths. Both here and for the remainder of the spectral code in
`pbrt`, wavelengths are specified in nanometers, which are of a magnitude
that gives easily human-readable values for the visible
wavelengths.

### 4.5.1 Spectrum Interface

We will find a variety of spectral representations useful in `pbrt`, ranging
from spectral sample values tabularized by wavelength to functional descriptions
such as the blackbody function. This brings us to our first interface
class, `Spectrum`. A `Spectrum` corresponds to a pointer
to a class that implements one such spectral representation.

`Spectrum` inherits from `TaggedPointer`, which handles the
details of runtime polymorphism. `TaggedPointer` requires that all
the types of `Spectrum` implementations be provided as template
parameters, which allows it to associate a unique integer identifier with
each type. (See Section B.4.4 for details of its
implementation.)

As with other classes that are based on `TaggedPointer`,
`Spectrum` defines an interface that must be implemented by all the
spectral representations. Typical practice in C++ would be for such an
interface to be specified by pure virtual methods in `Spectrum` and
for `Spectrum` implementations to inherit from `Spectrum` and
implement those methods. With the `TaggedPointer` approach, the
interface is specified implicitly: for each method in the interface, there
is a method in `Spectrum` that dispatches calls to the appropriate
type’s implementation. We will discuss the details of how this works for a
single method here but will omit them for other `Spectrum` methods and
for other interface classes since they all follow the same boilerplate.

The most important method that `Spectrum` defines is
`operator()`, which takes a single wavelength and returns
the value of the spectral distribution for that wavelength.

The corresponding method implementation is brief, though dense. A call to
`TaggedPointer::Dispatch()` begins the process of dispatching the
method call. The `TaggedPointer` class stores an integer tag along with the
object’s pointer that encodes its type; in turn, `Dispatch()` is able
to determine the specific type of the pointer at runtime. It then calls
the callback function provided to it with a pointer to the object, cast to
be a pointer to its actual type.

The lambda function that is called here, `op`, takes a pointer with
the `auto` type specifier for its parameter. In C++17, such a lambda
function acts as a templated function; a call to it with a concrete type
acts as an instantiation of a lambda that takes that type. Thus, the call
`(*ptr)(lambda)` in the lambda body ends up as a direct call to the
appropriate method.

`Spectrum` implementations must also provide a `MaxValue()`
method that returns a bound on the maximum value of the spectral
distribution over its wavelength range. This method’s main use in `pbrt` is for
computing bounds on the power emitted by light sources so that lights can
be sampled according to their expected contribution to illumination in the
scene.

### 4.5.2 General Spectral Distributions

With the `Spectrum` interface specified, we will start by defining a
few `Spectrum` class implementations that explicitly tabularize values
of the spectral distribution function. `ConstantSpectrum` is the
simplest: it represents a constant spectral distribution over all
wavelengths. The most common use of the `ConstantSpectrum` class in `pbrt` is
to define a zero-valued spectral distribution in cases where a particular
form of scattering is not present.

The `ConstantSpectrum` implementation is straightforward and we omit
its trivial `MaxValue()` method here. Note that it does not
inherit from `Spectrum`. This is another difference from using
traditional C++ abstract base classes with virtual functions—as far as
the C++ type system is concerned, there is no explicit connection between
`ConstantSpectrum` and `Spectrum`.

More expressive is `DenselySampledSpectrum`, which stores a spectral
distribution sampled at 1 nm intervals over a given range of integer
wavelengths .

Its constructor takes another `Spectrum` and evaluates that
spectral distribution at each wavelength in the range.
`DenselySampledSpectrum` can be useful if the provided spectral
distribution is computationally expensive to evaluate, as it allows
subsequent evaluations to be performed by reading a single value from
memory.

Finding the spectrum’s value for a given wavelength `lambda` is a
matter of returning zero for wavelengths outside of the valid range and
indexing into the stored values
otherwise.

While sampling a spectral distribution at 1 nm wavelengths gives sufficient
accuracy for most uses in rendering, doing so requires nearly 2 kB of
memory to store a distribution that covers the visible wavelengths.
`PiecewiseLinearSpectrum` offers another representation that is often
more compact; its distribution is specified by a set of pairs of values
where the spectral distribution is defined by
linearly interpolating between them; see Figure 4.17. For
spectra that are smooth in some regions and change rapidly in others, this
representation makes it possible to specify the distribution at a higher
rate in regions where its variation is greatest.

`PiecewiseLinearSpectrum`defines a spectral distribution using a set of sample values . A continuous distribution is then defined by linearly interpolating between them.

The `PiecewiseLinearSpectrum` constructor, not included here, checks
that the provided `lambda` values are sorted and then stores them and
the associated spectrum values in corresponding member variables.

Finding the value for a given wavelength requires first finding the pair of
values in the `lambdas` array that bracket it and then linearly interpolating
between them.

`PiecewiseLinearSpectrum`corner cases>> <<Find offset to largest

`lambdas`below

`lambda`and interpolate>> }

As with `DenselySampledSpectrum`, wavelengths outside of the specified
range are given a value of zero.

`PiecewiseLinearSpectrum`corner cases>>=

If `lambda` is in range, then `FindInterval()` gives the offset to
the largest value of `lambdas` that is less than or equal to
`lambda`. In turn, `lambda`’s offset between that wavelength and
the next
gives the linear interpolation parameter to use with the stored
values.

`lambdas`below

`lambda`and interpolate>>=

The maximum value of the distribution is easily
found using `std::max_element()`, which performs a linear search.
This function is not currently called in any performance-sensitive parts of
`pbrt`; if it was, it would likely be worth caching this value to avoid
recomputing it.

Another useful `Spectrum` implementation, `BlackbodySpectrum`,
gives the spectral distribution of a blackbody emitter at a specified
temperature.

The temperature of the blackbody in Kelvin is the constructor’s only parameter.

Because the power emitted by a blackbody grows so quickly with temperature
(recall the Stefan–Boltzmann law,
Equation (4.19)), the `BlackbodySpectrum`
represents a normalized blackbody spectral distribution where the maximum
value at any wavelength is 1. Wien’s displacement law,
Equation (4.20), gives the wavelength in meters where
emitted radiance is at its maximum; we must convert this value to nm before
calling `Blackbody()` to find the corresponding radiance value.

The method that returns the value of the distribution at a wavelength then
returns the product of the value returned by `Blackbody()` and the
normalization factor.

### 4.5.3 Embedded Spectral Data

`pbrt`’s scene description format provides multiple ways to specify spectral
data, ranging from blackbody temperatures to arrays of -value
pairs to specify a piecewise-linear spectrum. For convenience, a variety
of useful spectral distributions are also embedded directly in the `pbrt` binary, including ones that describe the emission profiles of various types
of light source, the scattering properties of various conductors, and the
wavelength-dependent indices of refraction of various types of glass. See
the online `pbrt` file format documentation for a list of all of them.

The `GetNamedSpectrum()` function searches through these spectra and
returns a `Spectrum` corresponding to a given named spectrum if it is
available.

A number of important spectra are made available directly through
corresponding functions, all of which are in a `Spectra` namespace.
Among them are `Spectra::X()`,
`Spectra::Y()`, and
`Spectra::Z()`, which return the
color matching curves that are described in Section 4.6.1,
and `Spectra::D()`, which returns a `DenselySampledSpectrum`
representing the D illuminant at the given temperature.

### 4.5.4 Sampled Spectral Distributions

The attentive reader may have noticed that although `Spectrum`
makes it possible to evaluate spectral distribution functions, it does not
provide the ability to do very much computation with them other than
sampling their value at a specified wavelength. Yet, for example,
evaluating the integrand of the reflection equation,
(4.14), requires taking the product of two spectral
distributions, one for the BSDF and one for the incident radiance function.

Providing this functionality with the abstractions that have been
introduced so far would quickly become unwieldy. For example, while the
product of two `DenselySampledSpectrum`s could be faithfully
represented by another `DenselySampledSpectrum`, consider taking the
product of two `PiecewiseLinearSpectrum`s: the resulting function would
be piecewise- quadratic and subsequent products would only increase its
degree. Further, operations between `Spectrum` implementations of
different types would not only require a custom implementation for each
pair, but would require choosing a suitable `Spectrum` representation
for each result.

`pbrt` avoids this complexity by performing spectral calculations at a set
of discrete wavelengths as part of the Monte Carlo integration that is
already being performed for image synthesis. To understand how this works,
consider computing the (non-spectral) irradiance at some point with
surface normal over some range of wavelengths of interest,
. Using Equation
(4.7), which expresses irradiance in terms of
incident radiance, and Equation (4.5), which
expresses radiance in terms of spectral radiance, we have

where is the incident spectral radiance at wavelength .

Applying the standard Monte Carlo estimator and taking advantage of the fact that and are independent, we can see that estimates of can be computed by sampling directions from some distribution , wavelengths from some distribution , and then evaluating:

Thus, we only need to be able to evaluate the integrand at the specified
discrete wavelengths to estimate the irradiance.
More generally, we will see that it is possible to express all the
spectral quantities that `pbrt` outputs as integrals over wavelength.
For example, Section 4.6 shows that
when rendering an image represented using RGB colors, each pixel’s color
can be computed by integrating the spectral radiance arriving at a pixel
with functions that model red, green, and blue color response. `pbrt` therefore uses only discrete spectral samples for spectral computation.

So that we can proceed to the implementation of the classes related to
sampling spectra and performing computations with spectral samples, we will
define the constant that sets the number of spectral samples here.
(Section 4.6.5 will discuss in more detail the trade-offs
involved in choosing this value.) `pbrt` uses 4 wavelength
samples by default; this value can easily be changed, though doing so
requires recompiling the system.

#### SampledSpectrum

The `SampledSpectrum` class stores an array of `NSpectrumSamples`
values that represent values of the spectral distribution at discrete
wavelengths. It provides methods that allow a variety of mathematical
operations to be performed with them.

Its constructors include one that allows providing a single value for all
wavelengths and one that takes an appropriately sized `pstd::span` of
per-wavelength values.

The usual indexing operations are also provided for accessing and setting each wavelength’s value.

It is often useful to know if all the values in a
`SampledSpectrum` are zero. For example, if a surface has zero
reflectance, then the light transport routines can avoid the computational
cost of casting reflection rays that have contributions that would
eventually be multiplied by zeros. This
capability is provided through a type conversion operator to
`bool`.

All the standard arithmetic operations on `SampledSpectrum` objects
are provided; each operates component-wise on the stored values. The
implementation of `operator+=` is below. The others are analogous and
are therefore not included in the text.

`SafeDiv()` divides two sampled spectra, but generates zero for any
sample where the divisor is zero.

In addition to the basic arithmetic operations, `SampledSpectrum` also provides
`Lerp()`,
`Sqrt()`,
`Clamp()`,
`ClampZero()`,
`Pow()`,
`Exp()`, and
`FastExp()` functions that
operate (again, component-wise) on `SampledSpectrum` objects; some of these
operations are necessary for evaluating some of the reflection
models in Chapter 9 and for evaluating volume
scattering models in Chapter 14. Finally,
`MinComponentValue()`
and
`MaxComponentValue()`
return the minimum and maximum of all the values, and
`Average()` returns their
average. These methods are all straightforward and are therefore not
included in the text.

#### SampledWavelengths

A separate class, `SampledWavelengths`, stores the wavelengths for
which a `SampledSpectrum` stores samples. Thus, it is important
not only to keep careful track of the `SampledWavelengths` that are
represented by an individual `SampledSpectrum` but also to not perform
any operations that combine `SampledSpectrum`s that have samples at
different wavelengths.

`u`>> <<Initialize

`lambda`for remaining wavelengths>>

`up`for th wavelength sample>>

To be used in the context of Monte Carlo integration, the
wavelengths stored in `SampledWavelengths` must be sampled from some
probability distribution. Therefore, the class stores
the wavelengths themselves as well as each one’s probability density.

The easiest way to sample wavelengths is uniformly over a
given range. This approach is implemented in the `SampleUniform()`
method, which takes a single uniform sample `u` and a range of
wavelengths.

`u`>> <<Initialize

`lambda`for remaining wavelengths>>

It chooses the first wavelength uniformly within the range.

`u`>>=

The remaining wavelengths are chosen by taking uniform steps `delta`
starting from the first wavelength and wrapping around if
`lambda_max` is passed. The result is a set of stratified wavelength
samples that are generated using a single random number. One advantage of
sampling wavelengths in this way rather than using a separate uniform
sample for each one is that the value of `NSpectrumSamples` can be
changed without requiring the modification of code that calls
`SampleUniform()` to adjust the number of sample values that are
passed to this method.

`lambda`for remaining wavelengths>>=

The probability density for each sample is easily computed, since the sampling distribution is uniform.

Additional methods provide access to the individual wavelengths and to all
of their PDFs. PDF values are returned in the form of a
`SampledSpectrum`, which makes it easy to compute the value of associated
Monte Carlo estimators.

In some cases, different wavelengths of light may follow different paths
after a scattering event. The most common example is when light undergoes
dispersion and different wavelengths of light refract to different
directions. When this happens, it is no longer possible to track multiple
wavelengths of light with a single ray. For this case, `SampledWavelengths`
provides the capability of terminating all but one of the wavelengths;
subsequent computations can then consider the single surviving wavelength
exclusively.

The wavelength stored in `lambda[0]` is always the survivor: there is
no need to randomly select the surviving wavelength so long as each
`lambda` value was randomly sampled from the same distribution as is
the case with `SampleUniform()`, for example. Note
that this means that it would be incorrect for
`SampledWavelengths::SampleUniform()` to always place `lambda[0]`
in a first wavelength stratum between `lambda_min` and
`lambda_min+delta`, `lambda[1]` in the second, and so
forth.

Terminated wavelengths have their PDF values set to zero; code that
computes Monte Carlo estimates using `SampledWavelengths` must
therefore detect this case and ignore terminated wavelengths accordingly.
The surviving wavelength’s PDF is updated to account for the termination
event by multiplying it by the probability of a wavelength surviving
termination, `1 / NSpectrumSamples`. (This is similar to how applying
Russian roulette affects the Monte Carlo estimator—see
Section 2.2.4.)

`SecondaryTerminated()` indicates whether `TerminateSecondary()`
has already been called. Because path termination is the only thing that
causes zero-valued PDFs after the first wavelength, checking the PDF values
suffices for this test.

We will often have a `Spectrum` and a set of wavelengths for
which we would like to evaluate it. Therefore, we will add a method to
the `Spectrum` interface that provides a
`Sample()` method that takes a set of wavelengths, evaluates its
spectral distribution function at each one, and returns a
`SampledSpectrum`. This convenience method eliminates the need for an
explicit loop over wavelengths with individual calls to
`Spectrum::operator()` in this common case. The implementations of this
method are straightforward and not included here.

#### Discussion

Now that `SampledWavelengths` and `SampledSpectrum` have been
introduced, it is reasonable to ask the question: why are they separate
classes, rather than a single class that stores both wavelengths and their
sample values? Indeed, an advantage of such a design would be that it
would be possible to detect at runtime if an operation was performed with
two `SampledSpectrum` instances that stored values for different
wavelengths—such an operation is nonsensical and would signify a bug in
the system.

However, in practice many `SampledSpectrum` objects are created during
rendering, many as
temporary values in the course of evaluating expressions involving spectral
computation. It is therefore worthwhile to minimize the object’s size, if
only to avoid initialization and copying of additional data.
While the `pbrt`’s CPU-based integrators do not store many
`SampledSpectrum` values in memory at the same time, the GPU rendering path
stores a few million of them, giving further motivation to minimize their
size.

Our experience has been that bugs from mixing computations at different
wavelengths have been rare. With the way that computation is structured in
`pbrt`, wavelengths are generally sampled at the start of following a ray’s
path through the scene, and then the same wavelengths are used throughout
for all spectral calculations along the path. There ends up being little
opportunity for inadvertent mingling of sampled wavelengths in
`SampledSpectrum` instances. Indeed, in an earlier version of the
system, `SampledSpectrum` did carry along a `SampledWavelengths`
member variable in debug builds in order to be able to check for that case.
It was eliminated in the interests of simplicity after a few months’
existence without finding a bug.