Spectral rendering

This chapter defines the spectral quantities used by HdRestir and relates the mathematical estimator to the CPU and GPU implementations. HdRestir stores four radiance samples per path. Each component corresponds to a wavelength between 360 and 830 nanometres. The four values remain associated with those wavelengths through emission, BSDF evaluation, absorption, and path throughput. Conversion to linear RGB occurs only when the renderer writes an image value or computes an RGB-derived diagnostic.

Spectral distributions and RGB values

A spectral radiance distribution \(L(\lambda)\) gives radiance per unit wavelength. Its units include an inverse wavelength term; if wavelength is measured in nanometres, one can write the spectral part as \(\mathrm{W\,sr^{-1}\,m^{-2}\,nm^{-1}}\). Integrating over a wavelength interval removes that inverse-nanometre term and gives radiance for that band.

An RGB value has a different meaning. It contains three coordinates in a specified color space, usually after integrating a spectrum against three response functions and applying a matrix transform. The mapping from spectra to three coordinates is many-to-one. Two different spectra that produce the same tristimulus coordinates are called metamers. Consequently, an RGB triplet cannot be inverted to recover one physically unique spectrum.

Wavelength-dependent refraction, absorption, fluorescence, and thin-film interference depend on the spectral variable directly. An RGB renderer can approximate some of these effects with fitted three-channel models, but it cannot generally represent a narrow feature that falls between its basis functions. Spectral transport retains wavelength as an integration variable, although its accuracy still depends on the spectra supplied for lights and materials.

From a spectrum to a pixel

The CIE 1931 standard observer dataset defines the three color-matching functions \(\bar x(\lambda),\bar y(\lambda),\bar z(\lambda)\). For a spectrum \(L(\lambda)\), the tristimulus values are

\(k\) is a normalization convention. HdRestir uses \(106.856\), an approximation to the integral of \(\bar y\) over the visible interval. It then applies the standard XYZ-to-linear-sRGB matrix. These operations produce scene-linear RGB; tone mapping and transfer-function encoding are separate post-processing operations.

The implementation does not tabulate the official CIE dataset. It evaluates the analytic fits published by Wyman, Sloan, and Shirley. The fits are compact and convenient for CPU/GPU parity, but they approximate the tabulated functions. The estimator below is unbiased for the functions actually evaluated by the program; approximation error in those functions is a separate source of error from Monte Carlo variance.

The integral is continuous, but a renderer cannot evaluate every wavelength. Choose wavelengths \(\lambda_i\) from a density \(p_\lambda\). The Monte Carlo estimate is

HdRestir’s four-wavelength packet

The implementation draws one \(u\sim U[0,1)\). Let \(\Delta=830-360=470\) nm and \(S=4\). It builds

The four samples are evenly separated, and the whole pattern receives a random cyclic shift. This is the central construction in hero-wavelength spectral sampling. HdRestir carries all four wavelengths through the path.

Even spacing reduces the chance that one path samples four wavelengths from a single narrow region. The random shift prevents a permanent fixed-grid bias.

The construction is easier to inspect by moving the initial random shift. Every marker moves by the same amount; crossing 830 nm wraps that marker to the start of the half-open interval.

flowchart LR
  U["one random u"] --> W["4 wrapped wavelengths"]
  W --> P["same geometric path"]
  P --> T["4 spectral throughputs"]
  T --> XYZ["CIE XYZ reconstruction"]
  XYZ --> RGB["linear RGB → display"]

Wavelength packets in the C++ implementation

Spectral sampling is already present in the first project commit, e6fe6bb. The following excerpts are literal selections from its source/core/spectrum.h; omitted members are not needed for this derivation.

#define SPECTRUM_SAMPLES 4
#define LAMBDA_MIN 360.0f
#define LAMBDA_MAX 830.0f

class SampledWavelengths {
public:
    float lambda[SPECTRUM_SAMPLES];

    // ...

    // u is a uniform random number in [0, 1)
    static SampledWavelengths SampleUniform(float u) {
        SampledWavelengths swl;
        swl.lambda[0] = LAMBDA_MIN + u * (LAMBDA_MAX - LAMBDA_MIN);
        for (int i = 1; i < SPECTRUM_SAMPLES; ++i) {
            swl.lambda[i] = swl.lambda[0] + (i * (LAMBDA_MAX - LAMBDA_MIN) / SPECTRUM_SAMPLES);
            if (swl.lambda[i] > LAMBDA_MAX) {
                swl.lambda[i] -= (LAMBDA_MAX - LAMBDA_MIN);
            }
        }
        return swl;
    }
};

The loop is the wrapped construction written mathematically above. With \(\Delta=470\) nm, the expression added for each successive component is \(i\Delta/4\). One random number controls the cyclic shift, so the code consumes one RNG value rather than four. The packet components are stratified but correlated.

The use of > rather than >= deserves a precise qualification. In the continuous derivation the probability of sampling exactly an endpoint is zero. Floating-point random numbers form a discrete set, so an exact 830.0F can in principle be retained by the implementation. This is a boundary convention, not evidence that the four components are independent. A future change that requires a strict half-open interval should wrap on >= and be covered by a boundary test.

The values transported at those wavelengths are stored separately:

class SampledSpectrum {
public:
    float values[SPECTRUM_SAMPLES];

    // ...

    SampledSpectrum operator*(const SampledSpectrum& s) const {
        SampledSpectrum ret;
        for (int i = 0; i < SPECTRUM_SAMPLES; ++i) ret.values[i] = values[i] * s.values[i];
        return ret;
    }

    // ...
};

Component-wise multiplication is appropriate because every factor refers to the same four wavelengths. A SampledSpectrum is not meaningful without the corresponding SampledWavelengths; the implementation passes the wavelength packet alongside spectral values rather than storing it inside every value.

Spectrum reconstruction in SpectrumToRGB

The CPU reconstruction first accumulates approximate CIE XYZ values:

float weight = (LAMBDA_MAX - LAMBDA_MIN) / SPECTRUM_SAMPLES;
// ...
weight /= 106.856f;

for (int i = 0; i < SPECTRUM_SAMPLES; ++i) {
    float l = lambda.lambda[i];
    float val = s.values[i];
    X += val * xFit_1931(l);
    Y += val * yFit_1931(l);
    Z += val * zFit_1931(l);
}

X *= weight;
Y *= weight;
Z *= weight;

(kLambdaMax-kLambdaMin)/kSpectrumSamples is \(\Delta/S=117.5\), which combines the reciprocal uniform PDF with the packet average. Dividing by 106.856F applies the chosen luminance normalization. The result is then transformed to linear sRGB:

float r =  3.2404542f * X - 1.5371385f * Y - 0.4985314f * Z;
float g = -0.9692660f * X + 1.8760108f * Y + 0.0415560f * Z;
float b =  0.0556434f * X - 0.2040259f * Y + 1.0572252f * Z;

// ...

r /= 1.198185f;
g /= 0.950366f;
b /= 0.907827f;

return GfVec3f(std::max(0.0f, r), std::max(0.0f, g), std::max(0.0f, b));

The three divisors white-balance the analytic reconstruction so a flat equal-energy spectrum maps approximately to (1,1,1). The final clamp removes negative linear-sRGB components. Negative components can be valid intermediate coordinates for spectra outside the sRGB gamut, so this clamp is a gamut handling choice rather than part of the CIE integral.

It also changes the statistical statement: XYZ before the clamp follows the unbiased Monte Carlo derivation; applying max(0,x) is nonlinear, therefore the expected clamped RGB value need not equal the clamp of the exact expected RGB value. The clamp prevents negative framebuffer values but can introduce small bias near the gamut boundary.

Where spectra enter the renderer

Scene materials and lights commonly arrive as RGB values. HdRestir “uplifts” those values to a smooth spectrum evaluated at the four active wavelengths. Multiplication is then component-wise: emission, BSDF value, and path throughput remain spectral until reconstruction.

The uplift implemented by RGBToSpectrum separates the achromatic component and adds three Gaussian lobes:

float minRGB = std::min(rgb[0], std::min(rgb[1], rgb[2]));
float white = minRGB;

GfVec3f rem(rgb[0] - white, rgb[1] - white, rgb[2] - white);

// ...

const float R_peak = 615.0f;
const float G_peak = 540.0f;
const float B_peak = 460.0f;
const float sigma = 20.0f;

for (int i = 0; i < SPECTRUM_SAMPLES; ++i) {
    float l = lambda.lambda[i];

    float val = white;

    // ...

    val += rem[0] * Gaussian(l, R_peak, sigma) * 1.295f;
    val += rem[1] * Gaussian(l, G_peak, sigma) * 1.556f;
    val += rem[2] * Gaussian(l, B_peak, sigma) * 1.382f;

    s.values[i] = std::max(0.0f, val);
}

This is a renderer-specific basis, not a physical inversion of RGB. The peak positions, widths, and scale factors are chosen to provide a useful round-trip. They do not guarantee that the reconstructed spectrum is the reflectance or emission spectrum of the original object, nor do they enforce all constraints required of measured reflectance data.

The GPU implementation does not exist at this point in the history. It is introduced by commit 4648adc and is discussed in the GPU chapters, where its constants can be compared with this original CPU implementation.

See the first implementation in spectrum.h and spectrum.cpp. For a broader reference, read PBRT’s radiometry, spectra, and color chapter.

Terminology used in the following chapters

The next chapter adds the other integral in the renderer: light arriving from directions and paths through the scene.