Incorrect illuminant when integrating spectral reflectance data to sRGB

[linusmossberg]

I found a problem when trying to verify my own solution to this that (at least) affects materials that uses spectral data to describe reflectance/transmittance properties. The problem is the following function:

static RGBSpectrum FromSampled(const Float *lambda, const Float *v, int n) {
    ...
    Float xyz[3] = {0, 0, 0};
    for (int i = 0; i < nCIESamples; ++i) {
        Float val = InterpolateSpectrumSamples(lambda, v, n, CIE_lambda[i]);
        xyz[0] += val * CIE_X[i];
        xyz[1] += val * CIE_Y[i];
        xyz[2] += val * CIE_Z[i];
    }
    Float scale = Float(CIE_lambda[nCIESamples - 1] - CIE_lambda[0]) /
                  Float(CIE_Y_integral * nCIESamples);
    xyz[0] *= scale;
    xyz[1] *= scale;
    xyz[2] *= scale;
    return FromXYZ(xyz);
}

Which corresponds to the Riemann sum of:

X = ∫(S(λ)·x(λ)·dλ) / ∫(y(λ)·dλ)
Y = ∫(S(λ)·y(λ)·dλ) / ∫(y(λ)·dλ)
Z = ∫(S(λ)·z(λ)·dλ) / ∫(y(λ)·dλ)

The problem for the reflectance case is that the reflectance S(λ) has to attenuate some illuminant I(λ) for these values to make sense, and the normalization factor should also take this into account:

X = ∫(S(λ)·I(λ)·x(λ)·dλ) / ∫(y(λ)·I(λ)·dλ)
Y = ∫(S(λ)·I(λ)·y(λ)·dλ) / ∫(y(λ)·I(λ)·dλ)
Z = ∫(S(λ)·I(λ)·z(λ)·dλ) / ∫(y(λ)·I(λ)·dλ)

The code implicitly uses a constant equal energy illuminant I(λ)=C that is eliminated, which is fine except that the function used to then transform the resulting XYZ tristimulus values to sRGB assumes XYZ_D65:

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];
}

This shifts the RGB values to more reddish hues in places where SPD's are used to describe reflectance, for example when specifying eta and k for conductive materials. The following scene uses spectral data to specify eta for usual soda-lime glass:

LookAt 3 4 1.5
       0 0 0
       0 0 1
Camera "perspective" "float fov" 25

Sampler "halton" "integer pixelsamples" 128
Integrator "path"
Film "image" "string filename" "glass-reflectance.png"
     "integer xresolution" [400] "integer yresolution" [400]

WorldBegin

LightSource "infinite" "rgb L" [1 1 1]
LightSource "distant"  "point from" [ -30 40  100 ] "rgb L" [10 10 10]

AttributeBegin
  Material "metal"

    "spectrum k" [ 352.9 0, 781.9 0 ]

    # https://refractiveindex.info/?shelf=glass&book=soda-lime&page=Rubin-clear
    "spectrum eta" [ 352.9 1.5444188160462
                     395.8 1.5377943708226
                     438.7 1.53297644448
                     481.6 1.5293470967242
                     524.5 1.5265302298432
                     567.4 1.5242862988783
                     610.3 1.5224568649146
                     653.2 1.5209337417404
                     696.1 1.5196410153786
                     739.0 1.5185241450581
                     781.9 1.5175431287183 ]

  Shape "sphere" "float radius" 1
AttributeEnd

WorldEnd

(FrConductor generalizes to dielectrics for k=0, so this is supposed to represent the specular reflectance part of glass)

This can be fixed by either using the D65 illuminant as I(λ) when applying the CIE CMF's in the reflectance case, or by simply multiplying the XYZ values with the D65 white point before applying the transformation to sRGB:

static RGBSpectrum FromSampled(const Float *lambda, const Float *v, int n, 
                               SpectrumType type = SpectrumType::Illuminant) {
    ...
    Float xyz[3] = {0, 0, 0};
    for (int i = 0; i < nCIESamples; ++i) {
        Float val = InterpolateSpectrumSamples(lambda, v, n, CIE_lambda[i]);
        xyz[0] += val * CIE_X[i];
        xyz[1] += val * CIE_Y[i];
        xyz[2] += val * CIE_Z[i];
    }
    Float scale = Float(CIE_lambda[nCIESamples - 1] - CIE_lambda[0]) /
                  Float(CIE_Y_integral * nCIESamples);
    xyz[0] *= scale;
    xyz[1] *= scale;
    xyz[2] *= scale;
    if (type == SpectrumType::Reflectance) {
        xyz[0] *= 0.95047f;
        xyz[2] *= 1.08883f;
    }
    return FromXYZ(xyz);
}

Note that this requires an additional SpectrumType argument that differentiates between reflectance and illuminant spectrums.