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.

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.

A controlled RGB comparison still requires two transport implementations with the same camera, light sources, material parameters, sampling budget, exposure, and output transform. The material sweep belongs to the implementation chapter, where its GGX parameters can be interpreted against the corresponding code. A future dispersion experiment should use a transmissive material with a documented wavelength-dependent IOR and a non-dispersive control.

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.