Add SpectralFilm.

This writes spectral images using the layout proposed in "An OpenEXR Layout
for Sepctral Images" by Fichet et al., https://jcgt.org/published/0010/03/01/.

src/pbrt/base/film.h

@@ -18,10 +18,11 @@ namespace pbrt {
 class VisibleSurface;
 class RGBFilm;
 class GBufferFilm;
+class SpectralFilm;
 class PixelSensor;
 
 // Film Definition
-class Film : public TaggedPointer<RGBFilm, GBufferFilm> {
+class Film : public TaggedPointer<RGBFilm, GBufferFilm, SpectralFilm> {
   public:
     // Film Interface
     PBRT_CPU_GPU inline void AddSample(Point2i pFilm, SampledSpectrum L,

src/pbrt/film.cpp

@@ -29,6 +29,7 @@
 #include <pbrt/util/transform.h>
 
 #include <algorithm>
+#include <cstring>
 
 namespace pbrt {
 
@@ -835,6 +836,213 @@ GBufferFilm *GBufferFilm::Create(const ParameterDictionary &parameters,
                                          writeFP16, alloc);
 }
 
+// SpectralFilm Method Definitions
+SpectralFilm::SpectralFilm(FilmBaseParameters p, Float lambdaMin, Float lambdaMax,
+                           int nBuckets, const RGBColorSpace *colorSpace,
+                           Float maxComponentValue, bool writeFP16, Allocator alloc)
+    : FilmBase(p),
+      colorSpace(colorSpace),
+      lambdaMin(lambdaMin),
+      lambdaMax(lambdaMax),
+      nBuckets(nBuckets),
+      maxComponentValue(maxComponentValue),
+      writeFP16(writeFP16),
+      pixels(p.pixelBounds, alloc) {
+    // Compute _outputRGBFromSensorRGB_ matrix
+    outputRGBFromSensorRGB = colorSpace->RGBFromXYZ * sensor->XYZFromSensorRGB;
+
+    filterIntegral = filter.Integral();
+    CHECK(!pixelBounds.IsEmpty());
+    filmPixelMemory += pixelBounds.Area() * (sizeof(Pixel) + 3 * nBuckets * sizeof(double));
+
+    // Allocate memory for the pixel buffers in big arrays. Note that it's
+    // wasteful (but convenient) to be storing three pointers in each
+    // SpectralFilm::Pixel structure since the addresses could be computed
+    // based on the base pointers and pixel coordinates.
+    int nPixels = pixelBounds.Area();
+    double *bucketWeightBuffer = alloc.allocate_object<double>(2 * nBuckets * nPixels);
+    std::memset(bucketWeightBuffer, 0, 2 * nBuckets * nPixels * sizeof(double));
+    AtomicDouble *splatBuffer = alloc.allocate_object<AtomicDouble>(nBuckets * nPixels);
+    std::memset(splatBuffer, 0, nBuckets * nPixels * sizeof(double));
+
+    for (Point2i p : pixelBounds) {
+        Pixel &pixel = pixels[p];
+        pixel.bucketSums = bucketWeightBuffer;
+        bucketWeightBuffer += NSpectrumSamples;
+        pixel.weightSums = bucketWeightBuffer;
+        bucketWeightBuffer += NSpectrumSamples;
+        pixel.bucketSplats = splatBuffer;
+        splatBuffer += NSpectrumSamples;
+    }
+}
+
+RGB SpectralFilm::GetPixelRGB(Point2i p, Float splatScale) const {
+    // Note: this is effectively the same as RGBFilm::GetPixelRGB
+
+    const Pixel &pixel = pixels[p];
+    RGB rgb(pixel.rgbSum[0], pixel.rgbSum[1], pixel.rgbSum[2]);
+    // Normalize _rgb_ with weight sum
+    Float weightSum = pixel.rgbWeightSum;
+    if (weightSum != 0)
+        rgb /= weightSum;
+
+    // Add splat value at pixel
+    for (int c = 0; c < 3; ++c)
+        rgb[c] += splatScale * pixel.rgbSplat[c] / filterIntegral;
+
+    // Convert _rgb_ to output RGB color space
+    rgb = outputRGBFromSensorRGB * rgb;
+
+    return rgb;
+}
+
+void SpectralFilm::AddSplat(Point2f p, SampledSpectrum L, const SampledWavelengths &lambda) {
+    // This, too, is similar to RGBFilm::AddSplat(), with additions for
+    // spectra.
+
+    CHECK(!L.HasNaNs());
+
+    // Convert sample radiance to _PixelSensor_ RGB
+    RGB rgb = sensor->ToSensorRGB(L, lambda);
+
+    // Optionally clamp sensor RGB value
+    Float m = std::max({rgb.r, rgb.g, rgb.b});
+    if (m > maxComponentValue)
+        rgb *= maxComponentValue / m;
+
+    // Spectral clamping and normalization.
+    Float lm = L.MaxComponentValue();
+    if (lm > maxComponentValue)
+        L *= maxComponentValue / lm;
+    L = SafeDiv(L, lambda.PDF()) / NSpectrumSamples;
+
+    // Compute bounds of affected pixels for splat, _splatBounds_
+    Point2f pDiscrete = p + Vector2f(0.5, 0.5);
+    Vector2f radius = filter.Radius();
+    Bounds2i splatBounds(Point2i(Floor(pDiscrete - radius)),
+                         Point2i(Floor(pDiscrete + radius)) + Vector2i(1, 1));
+    splatBounds = Intersect(splatBounds, pixelBounds);
+
+    // Splat both RGB and spectral bucket contributions.
+    for (Point2i pi : splatBounds) {
+        // Evaluate filter at _pi_ and add splat contribution
+        Float wt = filter.Evaluate(Point2f(p - pi - Vector2f(0.5, 0.5)));
+        if (wt != 0) {
+            Pixel &pixel = pixels[pi];
+
+            for (int i = 0; i < 3; ++i)
+                pixel.rgbSplat[i].Add(wt * rgb[i]);
+
+            for (int i = 0; i < NSpectrumSamples; ++i) {
+                int b = LambdaToBucket(lambda[i]);
+                pixel.bucketSplats[b].Add(wt * L[i]);
+            }
+        }
+    }
+}
+
+void SpectralFilm::WriteImage(ImageMetadata metadata, Float splatScale) {
+    Image image = GetImage(&metadata, splatScale);
+    LOG_VERBOSE("Writing image %s with bounds %s", filename, pixelBounds);
+    image.Write(filename, metadata);
+}
+
+Image SpectralFilm::GetImage(ImageMetadata *metadata, Float splatScale) {
+    // Convert image to RGB and compute final pixel values
+    LOG_VERBOSE("Computing final weighted pixel values");
+    PixelFormat format = writeFP16 ? PixelFormat::Half : PixelFormat::Float;
+
+    std::vector<std::string> imageChannels{{"R", "G", "B"}};
+    for (int i = 0; i < nBuckets; ++i) {
+        // The OpenEXR spectral layout takes the bucket center (and then
+        // determines bucket widths based on the neighbor wavelengths).
+        std::string lambda = StringPrintf("%.3fnm",
+                                          Lerp((i + 0.5f) / nBuckets, lambdaMin, lambdaMax));
+        // Convert any '.' to ',' in the number since OpenEXR uses '.' for
+        // separating layers.
+        std::replace(lambda.begin(), lambda.end(), '.', ',');
+
+        imageChannels.push_back("S0." + lambda);
+    }
+    Image image(format, Point2i(pixelBounds.Diagonal()), imageChannels);
+
+    std::atomic<int> nClamped{0};
+    ParallelFor2D(pixelBounds, [&](Point2i p) {
+        Pixel &pixel = pixels[p];
+
+        RGB rgb = GetPixelRGB(p, splatScale);
+
+        // Clamp to max representable fp16 to avoid Infs
+        if (writeFP16) {
+            for (int c = 0; c < 3; ++c) {
+                if (rgb[c] > 65504) {
+                    rgb[c] = 65504;
+                    ++nClamped;
+                }
+            }
+        }
+
+        Point2i pOffset(p.x - pixelBounds.pMin.x, p.y - pixelBounds.pMin.y);
+        image.SetChannels(pOffset, {rgb[0], rgb[1], rgb[2]});
+
+        // Set spectral channels. Hardcoded assuming that they come
+        // immediately after RGB, as is currently specified above.
+        for (int i = 0; i < nBuckets; ++i) {
+            Float c = 0;
+            if (pixel.weightSums[i] > 0) {
+                c = pixel.bucketSums[i] / pixel.weightSums[i] +
+                    splatScale * pixel.bucketSplats[i] / filterIntegral;
+                if (writeFP16 && c > 65504) {
+                    c = 65504;
+                    ++nClamped;
+                }
+            }
+            image.SetChannel(pOffset, 3 + i, c);
+        }
+    });
+
+    if (nClamped.load() > 0)
+        Warning("%d pixel values clamped to maximum fp16 value.", nClamped.load());
+
+    metadata->pixelBounds = pixelBounds;
+    metadata->fullResolution = fullResolution;
+    metadata->colorSpace = colorSpace;
+    metadata->strings["spectralLayoutVersion"] = "1.0";
+    // FIXME: if the RealisticCamera is being used, then we're actually
+    // storing "J.m^-2", but that isn't a supported value for
+    // "emissiveUnits" in the spec.
+    metadata->strings["emissiveUnits"] = "W.m^-2.sr^-1";
+
+    return image;
+}
+
+std::string SpectralFilm::ToString() const {
+    return StringPrintf(
+        "[ SpectralFilm %s lambdaMin: %f lambdaMax: %f nBuckets: %d writeFP16: %s maxComponentValue: %f ]",
+        BaseToString(), lambdaMin, lambdaMax, nBuckets, writeFP16, maxComponentValue);
+}
+
+SpectralFilm *SpectralFilm::Create(const ParameterDictionary &parameters, Float exposureTime,
+                                   Filter filter, const RGBColorSpace *colorSpace,
+                                   const FileLoc *loc, Allocator alloc) {
+    PixelSensor *sensor =
+        PixelSensor::Create(parameters, colorSpace, exposureTime, loc, alloc);
+    FilmBaseParameters filmBaseParameters(parameters, filter, sensor, loc);
+    bool writeFP16 = parameters.GetOneBool("savefp16", true);
+
+    if (!HasExtension(filmBaseParameters.filename, "exr"))
+        ErrorExit(loc, "%s: EXR is the only output format supported by the SpectralFilm.",
+                  filmBaseParameters.filename);
+
+    int nBuckets = parameters.GetOneInt("nbuckets", 16);
+    Float lambdaMin = parameters.GetOneFloat("lambdamin", Lambda_min);
+    Float lambdaMax = parameters.GetOneFloat("lambdamax", Lambda_max);
+    Float maxComponentValue = parameters.GetOneFloat("maxcomponentvalue", Infinity);
+
+    return alloc.new_object<SpectralFilm>(filmBaseParameters, lambdaMin, lambdaMax, nBuckets,
+                                          colorSpace, maxComponentValue, writeFP16, alloc);
+}
+
 Film Film::Create(const std::string &name, const ParameterDictionary &parameters,
                   Float exposureTime, const CameraTransform &cameraTransform,
                   Filter filter, const FileLoc *loc, Allocator alloc) {
@@ -845,6 +1053,9 @@ Film Film::Create(const std::string &name, const ParameterDictionary &parameters
     else if (name == "gbuffer")
         film = GBufferFilm::Create(parameters, exposureTime, cameraTransform, filter,
                                    parameters.ColorSpace(), loc, alloc);
+    else if (name == "spectral")
+        film = SpectralFilm::Create(parameters, exposureTime, filter,
+                                    parameters.ColorSpace(), loc, alloc);
     else
         ErrorExit(loc, "%s: film type unknown.", name);
 

src/pbrt/film.h

@@ -393,6 +393,127 @@ class GBufferFilm : public FilmBase {
     SquareMatrix<3> outputRGBFromSensorRGB;
 };
 
+// SpectralFilm Definition
+class SpectralFilm : public FilmBase {
+  public:
+    // SpectralFilm Public Methods
+    PBRT_CPU_GPU
+    bool UsesVisibleSurface() const { return false; }
+
+    PBRT_CPU_GPU
+    SampledWavelengths SampleWavelengths(Float u) const {
+        return SampledWavelengths::SampleUniform(u, lambdaMin, lambdaMax);
+    }
+
+    PBRT_CPU_GPU
+    void AddSample(Point2i pFilm, SampledSpectrum L, const SampledWavelengths &lambda,
+                   const VisibleSurface *, Float weight) {
+        // Start by doing more or less what RGBFilm::AddSample() does so
+        // that we can maintain accurate RGB values.
+
+        // Convert sample radiance to _PixelSensor_ RGB
+        RGB rgb = sensor->ToSensorRGB(L, lambda);
+
+        // Optionally clamp sensor RGB value
+        Float m = std::max({rgb.r, rgb.g, rgb.b});
+        if (m > maxComponentValue)
+            rgb *= maxComponentValue / m;
+
+        DCHECK(InsideExclusive(pFilm, pixelBounds));
+        // Update RGB fields in Pixel structure.
+        Pixel &pixel = pixels[pFilm];
+        for (int c = 0; c < 3; ++c)
+            pixel.rgbSum[c] += weight * rgb[c];
+        pixel.rgbWeightSum += weight;
+
+        // Spectral processing starts here.
+        // Optionally clamp spectral value. (TODO: for spectral should we
+        // just clamp channels individually?)
+        Float lm = L.MaxComponentValue();
+        if (lm > maxComponentValue)
+            L *= maxComponentValue / lm;
+
+        // The CIE_Y_integral factor effectively cancels out the effect of
+        // the conversion of light sources to use photometric units for
+        // specification.  We then do *not* divide by the PDF in |lambda|
+        // but take advantage of the fact that we know that it is uniform
+        // in SampleWavelengths(), the fact that the buckets all have the
+        // same extend, and can then just average radiance in buckets
+        // below.
+        L *= weight * CIE_Y_integral;
+
+        // Accumulate contributions in spectral buckets.
+        for (int i = 0; i < NSpectrumSamples; ++i) {
+            int b = LambdaToBucket(lambda[i]);
+            pixel.bucketSums[b] += L[i];
+            pixel.weightSums[b] += weight;
+        }
+    }
+
+    PBRT_CPU_GPU
+    RGB GetPixelRGB(Point2i p, Float splatScale = 1) const;
+
+    SpectralFilm(FilmBaseParameters p, Float lambdaMin, Float lambdaMax,
+                 int nBuckets, const RGBColorSpace *colorSpace,
+                 Float maxComponentValue = Infinity, bool writeFP16 = true,
+                 Allocator alloc = {});
+
+    static SpectralFilm *Create(const ParameterDictionary &parameters, Float exposureTime,
+                                Filter filter, const RGBColorSpace *colorSpace,
+                                const FileLoc *loc, Allocator alloc);
+
+    PBRT_CPU_GPU
+    void AddSplat(Point2f p, SampledSpectrum v, const SampledWavelengths &lambda);
+
+    void WriteImage(ImageMetadata metadata, Float splatScale = 1);
+
+    // Returns an image with both RGB and spectral components, following
+    // the layout proposed in "An OpenEXR Layout for Sepctral Images" by
+    // Fichet et al., https://jcgt.org/published/0010/03/01/.
+    Image GetImage(ImageMetadata *metadata, Float splatScale = 1);
+
+    std::string ToString() const;
+
+    PBRT_CPU_GPU
+    RGB ToOutputRGB(SampledSpectrum L, const SampledWavelengths &lambda) const {
+        LOG_FATAL("ToOutputRGB() is unimplemented. But that's ok since it's only used "
+                  "in the SPPM integrator, which is inherently very much based on "
+                  "RGB output.");
+        return {};
+    }
+
+  private:
+    PBRT_CPU_GPU
+    int LambdaToBucket(Float lambda) const {
+        DCHECK_RARE(1e6f, lambda < lambdaMin || lambda > lambdaMax);
+        int bucket = nBuckets * (lambda - lambdaMin) / (lambdaMax - lambdaMin);
+        return Clamp(bucket, 0, nBuckets - 1);
+    }
+
+    // SpectralFilm::Pixel Definition
+    struct Pixel {
+        Pixel() = default;
+        // Continue to store RGB, both to include in the final image as
+        // well as for previews during rendering.
+        double rgbSum[3] = {0., 0., 0.};
+        double rgbWeightSum = 0.;
+        AtomicDouble rgbSplat[3];
+        // The following will all have nBuckets entries.
+        double *bucketSums, *weightSums;
+        AtomicDouble *bucketSplats;
+    };
+
+    // SpectralFilm Private Members
+    const RGBColorSpace *colorSpace;
+    Float lambdaMin, lambdaMax;
+    int nBuckets;
+    Float maxComponentValue;
+    bool writeFP16;
+    Float filterIntegral;
+    Array2D<Pixel> pixels;
+    SquareMatrix<3> outputRGBFromSensorRGB;
+};
+
 PBRT_CPU_GPU
 inline SampledWavelengths Film::SampleWavelengths(Float u) const {
     auto sample = [&](auto ptr) { return ptr->SampleWavelengths(u); };