Computer Graphics 2022
ETH Zurich
Matias Turkulainen

Introduction

This website serves as a demo page for my work in the Computer Graphics course at ETH Zurich. I extended the minimalistic ray tracer Nori written in C++ with additional features to render the images shown here.

Motivation

My main motivation for the project was to tackle some of the challenging effects seen in images like the one below. I wanted to simulate volumetric scattering of light in a medium, difficult environmental maps, as well as realistic cameras with the thin lens equation.

Simple Features

Lens distortion

Lens distortion occurs in image formation with real lenses that deviate from the naive pinhole approximation. Light entering the optical lens off-center is distorted by the curvature of the lens. Radial distortion refers to the unequal bending of light depending on where light rays hit the lens relative to its optical center. Light further away from the optical center is distorted more. Radial distortion can be modelled using a polynomial radial distortion model as explained here (OpenCV). Mapping from undistorted pixels is determined by the radial coefficients K1, K2, and K3. Barell distortion occurs with positive coefficients and Pinchusion distortion occurs with negative coefficients.

 K1 = 0, K2 = 0, K3 = 0

 K1 = 1, K2 = 1, K3 = 0

 K1 = -1, K2 = -1, K3 = 0

/*Code implementation of polynomial distortion in perspective.cpp camera class
*/
float K1 = m_K1;
float K2 = m_K2;
float K3 = m_K3;

if (K1 != 0 || K2 != 0 || K3 != 0){
    // Original radius 
    float r = Vector2f(nearP.x() / nearP.z(), nearP.y() / nearP.z()).norm();
    // Polynomial Radial distortion
    nearP.x() *= (1 + K1 * pow(r, 2) + K2 * pow(r,4) + K3 * pow(r, 6));
    nearP.y() *= (1 + K1 * pow(r, 2) + K2 * pow(r,4) + K3 * pow(r, 6));
}

Depth of Field

Depth of field is a light transport phenomena that occurs through thin lenses. The size of the lens aperature and the focal length of the lens has an effect on where light is in focus on the image plane. Varying the aperature and the focal length changes what objects are in focus in the image.

 Varying lens radius with focal distance of 50. 100 Samples Per Pixel.

 
                                /* Code implementation of DoF in perspective.cpp camera class
*/
ray.o = m_cameraToWorld * Point3f(0, 0, 0);
ray.d = m_cameraToWorld * d;
ray.mint = m_nearClip * invZ;
ray.maxt = m_farClip * invZ;
ray.update();

// Implement Depth of Field
if(m_lensr > 0){ // Thin lens assumption
    //1. sample point on thin lens and scale by lensr
    Point2f pLens = m_lensr * Warp::squareToUniformDisk(apertureSample);

    // 2. compute intersection point on plane of focus. found by t*ray.z() = focal length
    float ft = m_focald/d.z();
    Point3f pFocus = d * ft;

    // 3. update ray for effect of lens. Ray origin and direction.  
    Point3f o(pLens.x(), pLens.y(), 0);
    ray.o = m_cameraToWorld * o;
    d = (pFocus - o).normalized();
    ray.d = m_cameraToWorld * d;
}
return Color3f(1.0f);

Image Textures

Texturing meshes in a scene is vital for any modern renderer. One way to texture meshes is to directly use images and to define a function that maps pixel colors to mesh locations. The method of image textures map pixel coordinates and associated RGB values in .png/.jpg images to object level uv-coordinates defined for a mesh in a scene. When a light ray intersects an object in the scene, the intersection point record stores a uv-coordinate value for the mesh intersection. This uv-coordinate is then mapped to a color value in a RGB image which is given to the shader such that an appropriate color can be rendered. Each mesh type in a renderer must define a unique mapping from 3D intersection points to uv-coordinates which can be used for texture mapping. The overall pipeline from light ray intersections to RGB values is illustrated below.

 f(x,y,z) -> (u,v); T(u,v) -> (R,G,B)

In Nori, I implement the image texture mapping fuction T(u,v). I map the uv-coordinate domain [0,1]² to pixel values on [height,width] with associated RGB color. I include the stb library for efficient image loading. I do not implement filtering, scale, or gamma correction, I do not cache or preload image textures in a scene (important for large scenes), and I only provide support for RGB Spectrum images.

I use high detail 4K texture maps of planets from the following source for my texture maps. See the code below for the full implementation and a few renders of image texture mapping of planets.

 

    /*
    Image Textures
    by: Matias Turkulainen
    */

    #include <nori/stb_image.h>
    #include <nori/texture.h>
    #include <filesystem/resolver.h>
    
    NORI_NAMESPACE_BEGIN
    
    class Image_texture : public Texture<Color3f> {
    public:
    
        Image_texture(const PropertyList &props) {
            m_path= props.getString("image_path", "");
            filesystem::path image_path = getFileResolver()->resolve(m_path);
            int width, height, channels;
            unsigned char *image = stbi_load(image_path.str().c_str(), &width, &height, &channels, 3);
            if(image == NULL) {
              printf("Error in loading the image\n");
              exit(1);
            }
    
            m_image = image;
            m_width = width;
            m_height = height;
            m_channels = channels;
        }
    
        // Map (u,v) to (r,g,b)
        Color3f eval(const Point2f & uv) {
            // map (u,v) [0,1] -> image pixels (x,y) [0,width/height]
    
            int x = int(uv.x() * m_width) ;
            int y = int((uv.y()) * m_height) ;
            
            float r = m_image[(y*m_width +x) * m_channels +0] /255.0 ;
            float g = m_image[(y*m_width +x) * m_channels +1] /255.0 ;
            float b = m_image[(y*m_width +x) * m_channels +2] /255.0;
            
            return Color3f(r, g, b);
        }
    
        std::string toString() const {
            return tfm::format(
                    "Image_texture[\n"
                    "  m_name = %s\n"
                    "]",
                    m_path // TODO: add rest of params toString
                );
        }
    
    protected:
        std::string m_path;
        unsigned char * m_image;
        unsigned m_width;
        unsigned m_height;
        unsigned m_channels;
    };
    
    NORI_REGISTER_CLASS(Image_texture, "image_texture")
    NORI_NAMESPACE_END

For validation, I attempted to recreate the same scene in Mitsuba3; however, I struggled with getting the correct lighting condition and camera transformation as in Nori. There is a rotation offset in the reference image.

Procedural Textures

A limitation with texture mapping is the need to store high resolution images which can still become blurry for close-up renders. Procedural textures attempt to combat this by defining the uv-coordinate to texture mapping with an explicit algorithm. Ideally, this allows for infinite resolution in your texture map, as you can query your procedural texture with the continuous uv-domain [0,1]. Finding a suitable algorithm for a desired texture is the crux of the challenge. Various algorithms exist for tessallation, such as the canonical checkerboard pattern (which I also implemented in Nori), and cellular texturing. More advanced mathematical functions, like the mandelbrot set, can also be defined as a procedural texture.

Another approach for generating textures algorithmically is to use a random texture generating function. The first random texture generator, referred to as the Image Synthesizer in the Ken Perlin's 1983 SIGGRAPH paper, was the Perlin noise function. It is a useful tool to naturally texture many objects in a scene without having to explicitly define image texture maps. Perlin textures have also been used to create virtual landscapes (publication date Feb, 2022) with varying parameters, e.g. height of hills or other scene elements. A notable example is the random scene generation in Minecraft.

The main idea behind the algorithm is a random function, referred to as a noise function, that maps 2D uv-coordinates to points in multiple random waves of varying amplitude and frequency which are summed together to create a texture output value. The key insight is that the mapping is made with a pseudo-random number generator which ensures that identical 2D uv-coordinates are mapped to the same output texture of the random function. This is a seeded random function and changing the seed with a scale attribute generates different random functions and thus output textures. Therefore, by controlling a single scale variable, different uv-coordinate to texture mappings can be easily created. Further control is given by explicitly defining the amplitude and frequency of the consequtively added waves so that the user can control the overall apperance and distribution of textures generated. More details regarding the algorithm and the following implementation can be found here: reference. For a slightly different treatment of Perlin noise and for a better overview of noise functions in procedural texture generation, refer to PBR.

My implementation of Perlin noise is controlled by a scale factor, persistence value, number of octaves, as well as an optional scaling factor for RGB colour output as seen in the .xml scene definition snippet below.

<!-- Sphere -->
 <mesh type="sphere">
    <point name="center" value="4.25,2.7,0"/>
    <float name="radius" value="2"/>
    <bsdf type="diffuse">
        <texture type="Procedural_perlin_texture" name="albedo">
            <integer name= "scale" value = "10"/>
            <float name= "persistence" value = "0.5"/>
            <integer name= "octaves" value = "8"/>
            <vector name = "color_scale" value = "10,1,1"/>
        </texture>
    </bsdf>
</mesh>

Number of octaves refers to the number of wave functions added together to generate the final texture. The higher the number of octaves, the more high frequency waves are added. Each octave, or individual wave, is defined by a frequency and amplitude. The frequency and amplitude are given by the following formula:

frequency = std::pow(2,i); amplitude = std::pow(persistence,i);

Where i is the octave number and persistence is a term that modulates the amplitude at each octave. Higher octaves have higher frequencies but lower amplitudes. Higher frequencies help in rendering fine details, like boundaries and edges in the generated texture patterns. The scale factor controls the pseudo-random number generator directly. I map uv-coordinates from the [0,1] domain to the interval [0, scale_factor]. These coordinates, modulated by the frequency at each octave, are the input to the noise function which deterministically generates random numbers on the interval [-1,1].

 // Noise function
float Noise(const float x, const float y){
int n;
n = x + y * 57;
n = (n << 13) ^ n;
float value = (1.0 - ((n * ((n * n * 15731) + 789221) + 1376312589) & 0x7fffffff) / 1073741824.0);
return (1.0 - ((n * ((n * n * 15731) + 789221) + 1376312589) & 0x7fffffff) / 1073741824.0);
}

The output of the random number generator is then used to interpolate a color value after summing the contribution of each octave wave. I interpolate with linear and cosine interpolation to create smooth texture transitions in the output color. The entire code and further implementation details can be seen below.

 
                /*
Procedural Perlin Noise Textures
by: Matias Turkulainen
*/

#include <nori/texture.h>
#include <filesystem/resolver.h>

NORI_NAMESPACE_BEGIN

class Procedural_perlin_texture : public Texture<Color3f> {
public:

    Procedural_perlin_texture(const PropertyList &props){
        m_persistence = props.getFloat("persistence", 0.5);
        m_octaves = props.getInteger("octaves", 8);
        m_scale = props.getInteger("scale", 10);
        m_color_scale = props.getVector3("color_scale", Vector3f(1));
    }
    
    Color3f eval(const Point2f & uv) {
        float total = 0.0;
        float p = m_persistence;
        int n = m_octaves - 1;

        int frequency;
        float amplitude;
        for (int i = 0; i < n; i++)
        {
            frequency = std::pow(2,i);
            amplitude = std::pow(p,i);
            total += amplitude * InterpolatedNoise(uv.x() * frequency * m_scale, uv.y() * frequency * m_scale);
            if (total >= 255){
                break;
            }
        }
        if (total >=0.0){
            return Color3f(total *m_color_scale.x(), total *m_color_scale.y(), total *m_color_scale.z());
        } else {return Color3f(0,0,0);}
    }

    float Noise(const float x, const float y){
        // Returns float in interval [-1,1]
        int n;
        n = x + y * 57;
        n = (n << 13) ^ n;
        float value = (1.0 - ((n * ((n * n * 15731) + 789221) +  1376312589) & 0x7fffffff) / 1073741824.0);
        return (1.0 - ((n * ((n * n * 15731) + 789221) +  1376312589) & 0x7fffffff) / 1073741824.0);
    }

    float SmoothedNoise(const float x, const float y){
        // Gradient smoothing
        float corners;
        float sides;
        float center;

        corners = ( Noise(x-1, y-1)+Noise(x+1, y-1)+Noise(x-1, y+1)+Noise(x+1, y+1) ) / 16;
        sides   = ( Noise(x-1, y)  +Noise(x+1, y)  +Noise(x, y-1)  +Noise(x, y+1) ) /  8;
        center  =  Noise(x, y) / 4;

        //std::cout<<"Smoothed noise is: "<<corners + sides + center<<std::endl;
        return corners + sides + center;
    }

    float InterpolatedNoise(const float x, const float y){
        // Grid computation
        int integer_X    = int(x);
        float fractional_X = x - integer_X;

        int integer_Y    = int(y);
        float fractional_Y = y - integer_Y;
        float v1 = SmoothedNoise(integer_X,     integer_Y);
        float v2 = SmoothedNoise(integer_X + 1, integer_Y);
        float v3 = SmoothedNoise(integer_X,     integer_Y + 1);
        float v4 = SmoothedNoise(integer_X + 1, integer_Y + 1);

        float i1 = Interpolate(v1 , v2 , fractional_X);
        float i2 = Interpolate(v3 , v4 , fractional_X);

        return Interpolate(i1 , i2 , fractional_Y);

    }

    float Interplate_linar (const float a, const float b, const float x){
        // Linear interpolation
        return  a*(1-x) + b*x;
    }

    float Interpolate(float a, float b, float x){
        // Cosine interpolation
        double ft = x * M_PI;
        double f = (1 - cos(ft)) * 0.5;
        return  a*(1-f) + b*f;
    }


    std::string toString() const {
        return tfm::format(
                "Procedural_perlin_texture[\n"
                "  m_name = %s\n"
                "]",
                m_persistence,
                m_octaves,
                m_scale
                );
    }

protected:
    int m_octaves;
    float m_persistence;
    int m_scale;
    Vector3f m_color_scale;
};

NORI_REGISTER_CLASS(Procedural_perlin_texture, "Procedural_perlin_texture")
NORI_NAMESPACE_END

 Scale = 10, Persistence = 0.5

 Scale = 100, Persistence = 0.5

I also explored interpolation techniques for smoother texture mapping. Linear interpolation lost details at edges and corners compared to higher order interpolation, like cosine interpolation.

 Linear vs Cosine interpolation, Scale = 50

Rendering on a High Performance Cluster

Some of the images seen already required over an hour to render on a M1 Macbook. To accelerate debugging and high resolution image creation, I attempted to compile Nori on a high performance cluster: Euler. To achieve this, some modifications were required. The major dependency on NanoGui was removed and a new main script was created that directly renders images bypassing previous GUI implementation steps. Redundant graphical dependencies were also removed from CMakeLists.txt file.

# Cluster Nori 
# by Matias Turkulainen
# This is the main CMake build file used to compile Nori Cluster Version
cmake_minimum_required (VERSION 2.8)
project("nori")

# This file included below takes care of compiling dependencies
# and setting compiler flags on different platforms
include("CMakeConfig.txt")

include_directories(ext)

# The following lines build the main executable. If you add a source
# code file to Nori, be sure to include it in this list.
add_executable(nori

    # Header files
    include/nori/cluster.h
    include/nori/bbox.h
    include/nori/bitmap.h
    include/nori/block.h
    include/nori/bsdf.h
    include/nori/bvh.h
    include/nori/camera.h
    include/nori/color.h
    include/nori/common.h
    include/nori/dpdf.h
    include/nori/frame.h
    include/nori/cluster.h
    include/nori/integrator.h
    include/nori/emitter.h
    include/nori/kdtree.h
    include/nori/mesh.h
    include/nori/object.h
    include/nori/parser.h
    include/nori/proplist.h
    include/nori/photon.h
    include/nori/ray.h
    include/nori/render.h
    include/nori/rfilter.h
    include/nori/sampler.h
    include/nori/scene.h
    include/nori/shape.h
    include/nori/texture.h
    include/nori/timer.h
    include/nori/transform.h
    include/nori/vector.h
    include/nori/warp.h

    # Source code files
    src/cluster.cpp
    src/main_cluster.cpp
    src/normals.cpp
    src/bitmap.cpp
    src/block.cpp
    src/bvh.cpp
    src/chi2test.cpp
    src/common.cpp
    src/consttexture.cpp
    src/checkerboard.cpp
    src/diffuse.cpp
    src/cluster.cpp
    src/independent.cpp
    src/mesh.cpp
    src/obj.cpp
    src/object.cpp
    src/parser.cpp
    src/perspective.cpp
    src/proplist.cpp
    src/render.cpp
    src/rfilter.cpp
    src/scene.cpp
    src/shape.cpp
    src/ttest.cpp
    src/warp.cpp
    src/microfacet.cpp
    src/photon.cpp
    src/mirror.cpp
    src/dielectric.cpp
    src/photonmapper.cpp
    src/sphere.cpp
    src/arealight.cpp
)

# The following lines build the warping test application
#add_executable(warptest
#  include/nori/warp.h
#  src/warp.cpp
#  src/warptest.cpp
#  src/microfacet.cpp
#  src/object.cpp
#  src/proplist.cpp
#  src/common.cpp
#)

#add_executable(tonemapper
#        include/nori/bitmap.h
#        src/bitmap.cpp
#        src/common.cpp
#        src/hdrToLdr.cpp)

# Nori depends on some libraries created in CMakeConfig.txt. The following two
# lines ensure that Nori is built *after* those libraries have been created.
add_dependencies(nori OpenEXR_p)
# add_dependencies(nori nanogui_p)
add_dependencies(nori tbb_p)
add_dependencies(nori pugixml)
#add_dependencies(warptest nori)
#add_dependencies(tonemapper nori)

# Link to dependency libraries
target_link_libraries(nori ${extra_libs})
#target_link_libraries(warptest ${extra_libs})
#target_link_libraries(tonemapper ${extra_libs})

# vim: set et ts=2 sw=2 ft=cmake nospell:

/*
Cluster version of Nori
by: Matias Turkulainen
*/

#include <nori/parser.h>
#include <nori/scene.h>
#include <nori/camera.h>
#include <nori/block.h>
#include <nori/timer.h>
#include <nori/bitmap.h>
#include <nori/sampler.h>
#include <nori/integrator.h>
#include <tbb/parallel_for.h>
#include <tbb/blocked_range.h>
#include <tbb/task_scheduler_init.h>
#include <filesystem/resolver.h>
#include <thread>

#include <nori/cluster.h>

using namespace nori;


static void renderBlock(const Scene *scene, Sampler *sampler, ImageBlock &block) {
    const Camera *camera = scene->getCamera();
    const Integrator *integrator = scene->getIntegrator();

    Point2i offset = block.getOffset();
    Vector2i size  = block.getSize();

    /* Clear the block contents */
    block.clear();

    /* For each pixel and pixel sample sample */
    for (int y=0; y<size.y(); ++y) {
        for (int x=0; x<size.x(); ++x) {
            for (uint32_t i=0; i<sampler->getSampleCount(); ++i) {
                Point2f pixelSample = Point2f((float) (x + offset.x()), (float) (y + offset.y())) + sampler->next2D();
                Point2f apertureSample = sampler->next2D();

                /* Sample a ray from the camera */
                Ray3f ray;
                Color3f value = camera->sampleRay(ray, pixelSample, apertureSample);

                /* Compute the incident radiance */
                value *= integrator->Li(scene, sampler, ray);

                /* Store in the image block */
                block.put(pixelSample, value);
            }
        }
    }
}

static void render(Scene *scene, const std::string &filename) {
    const Camera *camera = scene->getCamera();
    Vector2i outputSize = camera->getOutputSize();
    scene->getIntegrator()->preprocess(scene);

    /* Create a block generator (i.e. a work scheduler) */
    BlockGenerator blockGenerator(outputSize, NORI_BLOCK_SIZE);

    /* Allocate memory for the entire output image and clear it */
    ImageBlock result(outputSize, camera->getReconstructionFilter());
    result.clear();


    /* Do the following in parallel and asynchronously */
    std::thread render_thread([&] {
        tbb::task_scheduler_init init(-1);
        
        cout << "Rendering .. ";
        cout.flush();
        Timer timer;

        tbb::blocked_range<int> range(0, blockGenerator.getBlockCount());

        auto map = [&](const tbb::blocked_range<int> &range) {
            /* Allocate memory for a small image block to be rendered
                by the current thread */
            ImageBlock block(Vector2i(NORI_BLOCK_SIZE),
                camera->getReconstructionFilter());

            /* Create a clone of the sampler for the current thread */
            std::unique_ptr<Sampler> sampler(scene->getSampler()->clone());

            for (int i=range.begin(); i<range.end(); ++i) {
                /* Request an image block from the block generator */
                blockGenerator.next(block);

                /* Inform the sampler about the block to be rendered */
                sampler->prepare(block);

                /* Render all contained pixels */
                renderBlock(scene, sampler.get(), block);

                /* The image block has been processed. Now add it to
                    the "big" block that represents the entire image */
                result.put(block);
            }
        };

        /// Uncomment the following line for single threaded rendering
        // map(range);

        /// Default: parallel rendering
        tbb::parallel_for(range, map);

        cout << "done. (took " << timer.elapsedString() << ")" << endl;
    });


    /* Shut down the user interface */
    render_thread.join();

    /* Now turn the rendered image block into
        a properly normalized bitmap */
    std::unique_ptr<Bitmap> bitmap(result.toBitmap());

    /* Determine the filename of the output bitmap */
    std::string outputName = filename;
    size_t lastdot = outputName.find_last_of(".");
    if (lastdot != std::string::npos)
        outputName.erase(lastdot, std::string::npos);

    /* Save using the OpenEXR format */
    bitmap->saveEXR(outputName);

    /* Save tonemapped (sRGB) output using the PNG format */
    bitmap->savePNG(outputName);
}



int main(int argc, char **argv) {
    
    try {
        
        ImageBlock block(Vector2i(720, 720), nullptr);
        NoriCluster *cluster = new NoriCluster(block);

        // if file is passed as argument, handle it
        if (argc == 2) {
            std::string filename = argv[1];
            filesystem::path path(filename);

            if (path.extension() == "xml") {
                /* Render the XML scene file */
                    getFileResolver()->prepend(path.parent_path());

                std::unique_ptr<NoriObject> root(loadFromXML(argv[1]));

                /* When the XML root object is a scene, start rendering it .. */
                if (root->getClassType() == NoriObject::EScene)
                    render(static_cast<Scene *>(root.get()), argv[1]);

            } else {
                cerr << "Error: unknown file \"" << filename
                << "\", expected an extension of type .xml or .exr" << endl;
            }
        }

    } catch (const std::exception &e) {
        cerr << "Fatal error: " << e.what() << endl;
        return -1;
    }
    return 0;
}

Terminal output rendering on the ETH Euler Cluster (without making a batch job)

Advanced Features

... Whew, almost there ...

Environmental Map Emitter

Image based lighting refers to using a textured image as the envrionmental illumination source in a scene. This can be viewed as either an infinite area light that surrounds the entire scene or as a sphere that casts light into the scene from every direction. Images are favourable since high resolution images are able to capture high dynamic range effects and subtleties that are difficult to mimic with point lights, area emitters, or other light sources in a renderer. This approach is, unfortunately, very expensive since importance sampling of direct light sources is no longer effective since we effecively have an infinite area to sample. The light transport equation

must be approximated by Monte Carlo Sampling and a suitable importance sampling technique must be implemented which saves from integrating over the entire sphere for all scene intersection points. A sampling technique proportional to the outgoing illumination, or pixel wise intensity, on the environment map would help in lowering the variance on the incident radiance term. The main idea is to define a probability distribution on the pixel domain based on the environment map's luminance distribution (pixel wise intensities) and transforming this onto a probability distribution on a sphere. Then, the spherical distribution can be used to define a probability distribution on incident solid angles of scene intersection points. With this method, sampling incident solid angle directions proportional to the intensity of the environment map pixel intensities will greatly outperform and decrease the variance compared to naive uniform sampling of incident directions.

Below I walk through the implementation details for image based lighting in Nori. I closely follow the implementation illustrated in the paper Monte Carlo Rendering with Natural Illumination, especially their pseudocode for precomputation, 1D and 2D random sampling. I also took inspiration from the environment map implementation in the Mitsuba renderer.

Implementation Details

The first step in the algorithm is to precompute a piecewise-constant 2D distribution on the pixel domain proportional to the chosen environement map luminance or pixel wise intensitites f(x,y). The luminance is simply computed as the norm of the pixel RGB channels. The joint distribution of pixel luminance p(x,y) is then given by:

Sampling this distribution is not straightforward since the distributions p(x) and p(y) are not independent. Instead, we separate the joint distribution using the conditional property and marginalization. The marginal density p(y) can be expressed as:

This allows us to define the distribution of the remaining dimension as the following conditional distribution:

Therefore, the joint distribution p(x,y) simplifies to p(x,y) = p(y) x p(x|y) where both the marginal and conditional terms are precomputed at the environment map initialization. Sampling this 2D luminance distribution is then achieved by consecutive 1D piecewise sampling of the distribution p(y) followed by p(x|y). 1D sampling converts a random variable on the domain [0,1] (uv-domain) to the target pixel domain with the Inverse Method. However, our solution for the light transport equation requires a sampling distribution proportional to incident solid angle directions. It was only convenient to define the sampling distribution p(x,y) on the image, but this distribution must be mapped to solid angle directions. The following formula is used:

/*
Image Based Illumination: Environment Map Emitter
by: Matias Turkulainen

Image loading supports only OpenEXR format with the built in Nori/Bitmap class I/O

Reference code: 
1. https://libraopen.lib.virginia.edu/downloads/nv935289d
2. https://github.com/mitsuba-renderer/mitsuba/blob/master/src/emitters/envmap.cpp
*/

#include <nori/emitter.h>
#include <nori/bitmap.h>
#include <filesystem/resolver.h>
#include <nori/frame.h>
#include <iostream>
#include <fstream>
#include <cstdlib>

NORI_NAMESPACE_BEGIN

// Base Eigen matrix class used for bitmap 
typedef Eigen::Array<Color3f, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Base;

class Environment_map_emitter : public Emitter {
public:

    Environment_map_emitter(const PropertyList &props){
        m_path = props.getString("image_path", "");
        filesystem::path image_path = getFileResolver()->resolve(m_path);
        std::cout<<"Image path is: "<< m_path << std::endl;
        m_bitmap = Bitmap(image_path.str());
        
        m_rows = m_bitmap.rows();
        m_cols = m_bitmap.cols();

        // Convert bitmap RGB to intensity values
        Base m_intensity = compute_intensity_map();
        
        //-------------------------------------------//
        // Precomputation of conditional and marginal PDFs used for joint distribution of luminance p(x,y)

        // m_cdf is the conditional distribution p(x|y)
        m_cdf = Base(m_rows, m_cols + 1);
        // m_pdf 
        m_pdf = Base(m_rows, m_cols);
        
        // m_pdf_marginal is the marginal distribution p(x)
        m_pdf_marginal = Base(1, m_rows);
        m_cdf_marginal = Base(1, m_cols + 1);

        //precomputing the pdf and cdf values
        Base column_sum(1, m_rows); // rows = v
        for (int i = 0; i < m_pdf.rows(); ++i) {
            // conditional pdf and conditional cdf
            column_sum(0, i) = precompute(i, m_intensity, m_pdf, m_cdf); 
        }
        //marginal pdf and cdf
        precompute(0, column_sum, m_pdf_marginal, m_cdf_marginal);  
    }

    // luminance/intensity map f(x,y)
    Base compute_intensity_map(){
        Base intensity_map(m_rows, m_cols);
        for (int row = 0; row < m_rows; row++){
            for (int col = 0; col < m_cols; col++){
                intensity_map(row,col) = 
                                        std::pow(m_bitmap(row,col).r(),2) + 
                                        std::pow(m_bitmap(row,col).g(),2) + 
                                        std::pow(m_bitmap(row,col).b(),2); 
            }
        }
        return intensity_map;
    }

    // Precomputation function defined in Monte Carlo Rendering with Natural Illumination
    float precompute(int row, const Base &f, Base &pf, Base &Pf) const {
        float I = 0;
        int i;
        for(i= 0; i < f.cols(); i++)
            // Column sum
            I += f(row)(i);
        if (I == 0)
            return I;
        for(int j = 0; j < f.cols(); j++)
            //pdf definition
            pf(row, j) = f(row, j)/I;
        Pf(row, 0) = 0;
        for(i = 1; i < f.cols(); i++)
            // cumulative pdf definition
            Pf(row,i) = Pf(row, i-1) + pf(i-1);
        Pf(row,i) = 1;
        return I;
    }

    virtual Color3f eval(const EmitterQueryRecord &lRec) const {
        // return radiance (Color) at ref point
        // Pixel coordinates from incident angle direction vector
        Point2f coords = vectorToPixel(lRec.wi); 

        //TODO: Add maybe optional scaling on the output color
        return bilinearInterpolation(coords);
    }

    virtual Color3f sample(EmitterQueryRecord &lRec, const Point2f &sample) const {
        float u;
        float v;
        float u_pdf_marginal;
        float v_pdf_conditional;

        // Sample p(y) and then p(x|y) given y.
        sample1D(0, m_pdf_marginal, m_cdf_marginal, sample.x(), u, u_pdf_marginal);
        sample1D(u, m_pdf, m_cdf, sample.y(), v, v_pdf_conditional);

        // write samples to output file
        std::ofstream file;
        file.open("/Users/maturk/git/M1Nori/nori/scenes/project/environmentMap/sample/out_samples.txt", std::fstream::app);
        file << u<<","<<v<<"\n";
        file.close();

        Point2f pixel = Point2f(u, v);
        Vector3f w = pixelToVector(pixel);
        lRec.wi = w;
        lRec.shadowRay = Ray3f(lRec.ref, lRec.wi, Epsilon, 100000);

        // Jacobian transform from pixel to solid angle distribution
        float jacobian = (m_cols - 1) * (m_rows - 1) / (2 * std::pow(M_PI, 2) * Frame::sinTheta(lRec.wi));
        // pdf in solid angles
        float pdf_solid_angle = pdf(lRec) * jacobian;
        return eval(lRec) / pdf_solid_angle;
    }
    
    void sample1D(int rowNumber, const Base &pdf, const Base &cdf , const float &sample , float &x, float &prob) const {
        int i;
        for (i = 0; i < cdf.cols(); i++) {
            if ((sample >= cdf(rowNumber)(i)) && (sample < cdf(rowNumber)(i+1)))
                break;
        }

        float t = (cdf(rowNumber)(i+1) - sample) / (cdf(rowNumber)(i+1) - cdf(rowNumber)(i));
        x = (1-t) * i + t * (i+1);
        prob = pdf(rowNumber)(i);
    }

    virtual float pdf(const EmitterQueryRecord &lRec) const{
        // PDF p(x,y) is computed by p(x,y) = p_y(y) * p_x(x|y), i.e. marginal_y * conditional_x|y.
        // First convert lRec.wi vector to pixel value. 
        Point2f pixels = vectorToPixel(lRec.wi);
        int x = floor(pixels.x());
        int y = floor(pixels.y());

        return m_pdf_marginal(0)(y)*m_cdf(x)(y);
    }

    Point2f vectorToPixel(const Vector3f &vector) const {
        // Convert vector to spherical coordinates (theta, phi)
        Point2f coords = sphericalCoordinates(vector);
        float theta = coords.x();
        float phi = coords.y();
        // Convert spherical coordinates to uv-coordinates (u,v)
        // x: from [0, pi] to [0, m_rows -1]
        // y: from [-pi, pi] to [0, m_cols -1]
        float x = (m_rows - 1) * theta * INV_PI;
        float y = (m_cols - 1) * phi * 0.5 * INV_PI;
        
        return Point2f(x,y);
    }

    Vector3f pixelToVector(const Point2f &pixel) const{
        // from pixel (x,y) to spherical coordinates
        float theta = pixel.x() * M_PI /(m_rows - 1);
        float phi = pixel.y() * 2 * M_PI / (m_cols - 1);
        //spherical coordinates
        return Vector3f(sin(theta) * cos(phi), sin(theta)*sin(phi), cos(theta)).normalized();
    }



    Color3f bilinearInterpolation(Point2f coords) const {
        int x1 = floor(coords.x());
        int y1 = floor(coords.y());
        int x2 = x1 + 1;
        int y2 = y1 + 1;
        Color3f Q11 = 0.0;
        //TODO: do we need all the checks? investigate
        if (x1 >= 0 && x1 < m_rows && y1 >= 0 && y1 <m_cols)
            Q11 = m_bitmap(x1, y1); 
        Color3f Q12 = 0.f;
        if (x1 >= 0 && x1 < m_rows && y2 >= 0 && y2 < m_cols)
            Q12 = m_bitmap(x1, y2); 
        Color3f Q21 = 0.f;
        if (x2 >= 0 && x2 < m_rows && y1 >= 0 && y1 < m_cols)
            Q21 = m_bitmap(x2, y1); 
        Color3f Q22 = 0.f;
        if (x2 >= 0 && x2 < m_rows && y2 >= 0 && y2 < m_cols)
            Q22 = m_bitmap(x2, y2); 

        auto w11 = (x2-coords.x())*(y2-coords.y())/((x2-x1)*(y2-y1));
        auto w12 = (x2-coords.x())*(coords.y()-y1)/((x2-x1)*(y2-y1));
        auto w21 = (coords.x()-x1)*(y2-coords.y())/((x2-x1)*(y2-y1));
        auto w22 = (coords.x()-x1)*(coords.y()-y1)/((x2-x1)*(y2-y1));

        return w11*Q11 + w12*Q12 + w21*Q21 + w22*Q22;

    }

    virtual std::string toString() const {
        return "Environment_map_emitter[]";
    }

protected:
    std::string m_path;
    Bitmap m_bitmap;
    Bitmap m_intensity;
    int m_cols;
    int m_rows;
    Base m_pdf;
    Base m_cdf;
    Base m_pdf_marginal; 
    Base m_cdf_marginal;
    
};

NORI_REGISTER_CLASS(Environment_map_emitter, "Environment_map_emitter")
NORI_NAMESPACE_END

Evaluation

Here are some renders of my implementation compared to the Mitsuba 3 Physically Based Renderer:

Samples Per Pixel: 100
HDRI Image Reference:  https://polyhaven.com/a/satara_night

Samples Per Pixel: 100
HDRI Image Reference:  https://polyhaven.com/a/studio_garden

There are some minor issues with my implementation compared to the reference Mitsuba3 renders. I suspect it is due to my importance sampling in the pixel (x,y) domain not accounting for longitude-latitude distortion that occurs when mapping between different distribution domains.

Homogeneous Participating Media

Participating media in graphics rendering refers to scattering effects occuring between surfaces due to non empty space. This includes scattering effects due to the atmosphere, water, or various other mediums. Unlike BSDFs that characterize reflectance on single shape surfaces, scattering media are assumed to be so numerous that their scattering properties are represented by statistical distributions. Volume rendering is the extension of simple path integration, defined by the light transport equation in the previous section, in a renderer that enables light transport in participating media accounting for the various losses and additions in light radiance along a ray due to absorption and emission in a media. The goal of estimating the effect of participating media is to sample the following integral form of the light transport equation, referred to as the Equation of Transfer, which accounts for absorption, emission, and in-scattering of radiance along a ray path:

where L_i is the incident radiance at a point P in space which is equal to the reflected radiance back along the ray from a surface intersection point L_o and the added radiance along the ray due to volume scattering (emission, out-scattering and attenuation, and in-scattering) up to the surface intersection point. L_s refers to source radiance which accounts for the total added radiance per unit distance due to in-scattering and emission along a ray. Since most participating material are non-emissive, we assume that the source term only accounts for in-scattering, or the increased radiance due to scattering from other directions. The source term is given by

where σ_s is the in-scattering probability per unit distance, p is the phase function that describes the angular distribution of scattered radiation at a point (similar to BSDFs), and L_i is the incident radiance at that particular point in space integrated over the entire sphere of incident directions ω. Both terms in the Equation of Transport are attenuated by the beam transmittance function T_r which describes the fraction of radiance transmitted between two points.

Like most ray tracing methods, Monte Carlo sampling provides an unbiased estimate of the transport equation. Monte Carlo sampling the Equation of Transport results in two scenarios: 1. the scene intersection point radiance term L_o, occuring at a distance t away from the point p, is sampled and estimated and 2. a sample is drawn on the interval [0, t] and the participating media integral term is estimated. Here I only consider Homogeneous participating media which are constant density media bounded by surfaces, either within the entire scene bounding box or within shape primitives in a scene.

Implementation Details

The main implementation challenges are to properly define the sampling methods used for each interaction type and to develop a new type of renderer, referred to as a Volume Renderer, that casts rays with incremental step-sizes proportional to the density of the media. A participating media can be defined by a few key parameters, primarily by the absorption coefficient σ_a and the scattering coefficient σ_s which define the density of absorbed and scattered radiance per unit length respectively. These coefficients also define the extinction coefficient σ_t = σ_a + σ_s and the albedo α = σ_s / σ_t. The extinction coefficient, also referred to as the transmission coefficient, σ_t is used to define the beam transmittance between points in a media. For homogeneous material, the transmittance function has an easy closed form expression given by Beer's law:

where d is the distance between points p^' and p. Since the transmittance attenuates the contribution of radiance along ray paths due to the media, sampling step sizes proportional to the transmittance allows importance sampling of the medium. The free-flight distance, or the step-size, defines the distance to the next interaction point in a scene and I randomly draw steps proportional to the transmittance function of a homogeneous media. Step sizes t for the Beer's transmittance are given by:

where ξ is a random number in the range [0,1]. I extend a multiple-importance sampling (MIS) ray tracer implementation with incremental step sizes proportional to the transmittance distribution of the medium for each casted ray to create a volume renderer. Sampling surface interactions occurs when the step size t, defined above, exeeds the distance to the nearest scene intersection point. I also importance sample the phase function p(ω, p^', p) which defines the scattering directions ω_o within a media. For a homogeneous medium I assume scattering is isotropic such that the distribution of scattered radiance is given by uniform sphere sampling:

i.e. the distribution of scattered radiance is independent of incident direction (p, p_'). Russian roulette early stopping of ray casting is also implemented. Below are the class constuctors that are used to define a homogeneous medium in Nori.

medium.h
/* 
    Medium Class Header
    by Matias Turkulainen
*/

#if !defined(__NORI_MEDIUM_H)
#define __NORI_MEDIUM_H

#include <nori/object.h>
#include <nori/bbox.h>
#include <nori/sampler.h>
#include <nori/common.h>

NORI_NAMESPACE_BEGIN

struct MediumQueryRecord {

        // Medium query ref. Previous origin. 
        Point3f ref;

        // Medium intersection point
        Point3f p; 

        // Maximum distance i.e. scene intersection point
        float tMax;

        // out direction vector
        Vector3f wo;

        // Incident direction vector
        Vector3f wi;

        // PDF 
        float pdf;

        /// Create an unitialized query record
        MediumQueryRecord() : tMax(0), wo(0.0){ }

        MediumQueryRecord(const Point3f &ref) : p(ref) { }
        
};

class Medium : public NoriObject {

public:
    
    virtual ~Medium() { }

    virtual EClassType getClassType() const override { return EMedium; }

    // Get transmittance Tr 
    virtual float Tr(const Point3f &x, const Point3f &y) const = 0;

    // Sample medium
    virtual float sample(MediumQueryRecord &lRec, const Point2f &sample) const = 0;

    virtual Color3f eval(const MediumQueryRecord &lRec) const = 0;

    // Assume phase functions are sampled perfectly, set pdf to phaseFunction
    virtual float pdf(const MediumQueryRecord &lRec) const {
        return phaseFunction(lRec);
    }

    // Sample free flight distance 
    virtual float sampleTr(float sample) const = 0;

    // Sample phase functions
    // Isotropic
    virtual float phaseFunction(const MediumQueryRecord &lRec) const = 0;

    // Optional custom phase function depending on wo and wi
    float phaseFunction(const Vector3f &wo, const Vector3f &wi);

    // Henyey and Greenstein phase function
    inline float PhaseHG(float cosTheta, float g) {
        float denom = 1 + g * g + 2 * g * cosTheta;
        return INV_FOURPI * (1 - g * g) / (denom * std::sqrt(denom));
    }
    void setShape(Shape * shape) {
        m_shape = shape;
    }

protected:
    /// Pointer to the shape if the emitter is attached to a shape
    Shape * m_shape = nullptr;
    // Absorption coefficient    
    float m_sigma_a;

    //Scattering coefficient
    float m_sigma_s;

    //Extinction coefficient
    float m_sigma_t;

    //Albedo
    float m_albedo;
};


NORI_NAMESPACE_END

#endif /* __NORI_INTEGRATOR_H */

medium.cpp 
/*
Homogeneous Participating medium
by Matias Turkulainen
*/

#include <nori/medium.h>
#include <nori/warp.h>

NORI_NAMESPACE_BEGIN


class Homogeneous_participating_medium : public Medium {
public:

    Homogeneous_participating_medium(const PropertyList &props){
        m_sigma_a = props.getFloat("sigma_a", 0.1f);
        m_sigma_s = props.getFloat("sigma_s", 0.1f);
        m_sigma_t = m_sigma_a + m_sigma_s;
        m_albedo = m_sigma_s / m_sigma_t;

    }

    virtual ~Homogeneous_participating_medium() { }

    virtual Color3f eval(const MediumQueryRecord &lRec) const {
		return Color3f(1.0);
	}
    
    virtual float sample(MediumQueryRecord &lRec, const Point2f &sample) const {
		// Isotropic unifrom sampling results in unifrorm spherical sampling//
		lRec.wo = Warp::squareToUniformSphere(sample);
        // sample should give (Tr(t) * sigma_s * phase_function_value(w')) / (distance_sampling_pdf(t) * direction_sampling_pdf(w')) which evaluates to albedo for homogenous material
        if (m_albedo >= 0){return m_albedo;} 
        else {return 1.0;}
		
	}

    virtual float Tr(const Point3f &x, const Point3f &y) const {
        // Homogeneous transmission
        if (m_sigma_t > 0)
			return std::exp(-m_sigma_t * (y - x).norm());
		else
			return 1.f;
    }

    virtual float pdf(const MediumQueryRecord &lRec) const {
        return phaseFunction(lRec);
    }

    virtual float phaseFunction(const MediumQueryRecord &lRec) const {
		// Isotropic medium phase function
		return 1 / (4 * M_PI);
	}


    std::string toString() const {
		return tfm::format("Homogeneous_participating_medium[\n"
				"  m_sigma_a = %s,\n"
				"  m_sigma_s = %s,\n"
				"]", m_sigma_a, m_sigma_s);
	}


protected:
    float m_sigma_a;
    float m_sigma_s;
    float m_sigma_t;
	float m_albedo;

};

NORI_REGISTER_CLASS(Homogeneous_participating_medium, "Homogeneous_participating_medium")
NORI_NAMESPACE_END

Below is my implementation of a Volume Renderer with multiple importance sampling of medium transmittance, phase function, interaction type (medium or material sampling), as well as emitter sampling.

/*
Volumetric Path Integrator
by Matias Turkulainen
*/

#include <nori/integrator.h>
#include <nori/emitter.h>
#include <nori/scene.h>
#include <nori/bsdf.h>
#include <nori/sampler.h>
#include <nori/medium.h>
#include <nori/warp.h>
#include <nori/shape.h>

NORI_NAMESPACE_BEGIN

class Volumetric_path_integrator : public Integrator {
public:

    Volumetric_path_integrator(const PropertyList &props) {}
    
    Color3f Li(const Scene *scene, Sampler *sampler, const Ray3f &ray) const override {
        // Recursive Monte Carlo ray-tracing in a volume medium with russian roulette early stopping
        Color3f out_color(0.0);
        Color3f throughput = 1.0;
        float successProbability;
        Color3f brdf;
        float w_mats = 1.0;
        float w_ems = 1.0;
        
        Medium *medium = scene->getMedium()[0];

        Ray3f traceRay = ray; // Ray which is traced recursively in scene. Terminated with russian roulette. 
        Intersection its; 
        
        // Check if this ray hits the scene at all or if its intersection point is at infinity
        bool bounded = scene->rayIntersect(traceRay, its);

        // MIS + Russian roulette with medium
        while (true) { 
            MediumQueryRecord mi;
            float t_max;
            if (bounded)
                t_max = (traceRay.o - its.p).norm();
            else {t_max = 100000.0;} // its.p = infinity

            mi.tMax = t_max;

            // Sample free path distance, pdf is proportional to Tr
            float t = medium->sampleTr(sampler->next1D());
            
            // Sample medium
            if(t <= t_max - Epsilon){
                mi.p = traceRay.o + t * traceRay.d.normalized(); // move to medium location defined by t-step.
                float medium_color = medium->sample(mi, sampler->next2D()); // sample medium color at this point. sets wo to random sphere direction and medium_color = albedo/pdf(phaseFunction)
                
                float mats_pdf = medium->pdf(mi); // pdf at medium point. proportional to phaseFunction

                // Emitter Light Sampling
                // Sample lights at mi.p() location
                const Emitter *light = scene->getRandomEmitter(sampler->next1D());
                EmitterQueryRecord lRec(mi.p); // set reference point as point in medium
                Color3f Li = light->sample(lRec, sampler->next2D()) * scene->getLights().size(); // Sample light source, returns color weighted by num lights in scene
                float tr = medium->Tr(mi.p, lRec.p); // Transmittance between light source and point in medium
                float ems_pdf = light->pdf(lRec); // update emitter pdf

                if (!scene->rayIntersect(lRec.shadowRay)){ // check if occluded
                    throughput *= medium_color;
                    out_color += throughput * tr  * Li * mats_pdf; // Contribution due to light source
                    //std::cout<< "tr is "<<tr<<std::endl;
                    mi.tMax = lRec.shadowRay.maxt;
                }
                else{ // Light is occluded. Only account for medium color
                    throughput *= medium_color;
                }

                // Russian Roulett early stopping
                successProbability = fmin(throughput.x(), 0.99);
                if (sampler->next1D() > successProbability)
                    return out_color;
                else { // Continue to next ray
                    throughput /= successProbability;
                }
                traceRay = Ray3f(mi.p, mi.wo.normalized()); // New ray in the direction of phase function wo
                
                if(scene->rayIntersect(traceRay, its)){
                    if (its.mesh->isEmitter()) {
                        EmitterQueryRecord lRec = EmitterQueryRecord(traceRay.o, its.p, its.shFrame.n);
                        float pdf_ems = its.mesh->getEmitter()->pdf(lRec);
                        if (mats_pdf > 0 && pdf_ems + mats_pdf > 0){
                            w_mats = mats_pdf / (ems_pdf + mats_pdf);     
                        }
                        else {w_mats = mats_pdf;}
                        out_color += throughput * w_mats * its.mesh->getEmitter()->eval(lRec) * medium->Tr(traceRay.o, its.p);
                    }   
                }
            }

            // stop early if no intersection in the scene
            else if (!bounded) {
                return out_color;
            }

            else { // Surface intersection
                    // we are at point its.p which is a scene intersection point
                mi.p = traceRay.o + mi.tMax * traceRay.d.normalized();
                if (its.mesh->isEmitter())
                    {
                        EmitterQueryRecord lRec(traceRay.o, its.p, its.shFrame.n);
                        // use previous iteration w_mats calculation and add radiance mdulated by transmittance between intersection point and origin of traceRay
                        Color3f addition = throughput * w_mats * its.mesh->getEmitter()->eval(lRec) * medium->Tr(traceRay.o, its.p);
                        out_color += throughput * w_mats * its.mesh->getEmitter()->eval(lRec) * medium->Tr(traceRay.o, its.p);
                    }
        
                // Sample random emitter in scene
                EmitterQueryRecord lRec(its.p);
                const Emitter *light = scene->getRandomEmitter(sampler->next1D());
                lRec.ref = its.p; 
                Color3f Li = light->sample(lRec, sampler->next2D()) * scene->getLights().size(); // weighted irradiance 
                float ems_pdf = light->pdf(lRec); // update emitter pdf
        
                if (!scene->rayIntersect(lRec.shadowRay)){ // check if occluded
                    float theta = Frame::cosTheta(its.shFrame.toLocal(lRec.wi));
        
                    BSDFQueryRecord bRec(its.toLocal(-traceRay.d), its.shFrame.toLocal(lRec.wi), ESolidAngle);
                    float mats_pdf = its.mesh->getBSDF()->pdf(bRec); // update bsdf pdf
                    // Calcuate w_ems weight term:
                    if (ems_pdf > 0 && ems_pdf + mats_pdf > 0){
                        w_ems = ems_pdf / (ems_pdf + mats_pdf);   
                    }
                    else {w_ems = ems_pdf;}
                    //Update tMax
                    mi.tMax = lRec.shadowRay.maxt;
                    //Update ems_color according to weight, accounting for transmittance between o() and its.p()
                    out_color += throughput * medium->Tr(traceRay.o, lRec.p) * w_ems * its.mesh->getBSDF()->eval(bRec) * Li * theta;
                }
        
                // Russian Roulette early stopping
                successProbability = fmin(throughput.x(), 0.99);
                if (sampler->next1D() > successProbability)
                    return out_color;
                else { // Continue to next ray
                    throughput /= successProbability;
                }
        
                // query bsdf at scene intersection point
                BSDFQueryRecord bRec(its.shFrame.toLocal(-traceRay.d)); 
                bRec.uv = its.uv;
                
                brdf = its.mesh->getBSDF()->sample(bRec, sampler->next2D());
                throughput *= brdf; //(wi dot n)*bsdf/pdf
        
                // Store old mats_pdf
                float mats_pdf =  its.mesh->getBSDF()->pdf(bRec);
        
                traceRay = Ray3f(its.p, its.toWorld(bRec.wo)); // update traceRay,

                bool bounded = scene->rayIntersect(traceRay, its);

                if (!scene->rayIntersect(traceRay, its))
                    return out_color;
                
                if (its.mesh-> isEmitter()) {
                    EmitterQueryRecord lRec = EmitterQueryRecord(traceRay.o, its.p, its.shFrame.n);
                    float ems_pdf = its.mesh->getEmitter()->pdf(lRec);
        
                    if (mats_pdf > 0 && ems_pdf + mats_pdf > 0){
                        w_mats = mats_pdf / (ems_pdf + mats_pdf);     
                    }
                    else {w_mats = mats_pdf;}
                }
        
                // BUG! Needed to check measure values!!!
                if (bRec.measure == EDiscrete) {
                    // Lastly check for delta distribution point lights 
                    w_mats = 1.0;
                }
            }
        }
        return out_color;
    }

    std::string toString() const {
        return "Volumetric_path_integrator[]";
    }
};

NORI_REGISTER_CLASS(Volumetric_path_integrator, "Volumetric_path_integrator");
NORI_NAMESPACE_END

Evaluation

Here are some renders of my implementation compared to the Mitsuba 3 Physically Based Renderer:

Samples Per Pixel: 256
<medium type="Homogeneous_participating_medium">
<float name="sigma_a" value="0.001"/>
<float name="sigma_s" value="0.001"/>
</medium>

Samples Per Pixel: 256
Top Row
<medium type="Homogeneous_participating_medium">
<float name="sigma_a" value="0.1"/>
<float name="sigma_s" value="0.2"/>
</medium>
Bottom Row
<medium type="Homogeneous_participating_medium">
<float name="sigma_a" value="0.1"/>
<float name="sigma_s" value="0.5"/>
</medium>

Final Render

Here I attempt to recreate the motivational image. The sea floor is simulated with colored random perlin textures. I use my volume renderer to simulate absorption and scattering inside water. I simulate radial distortion and the depth of field effect with the underwater perspective camera. I also use an HDRI environment map of a sea to give more texture to the background.

Acknowledgements

Thanks to all the teachers assistants and people at Disney Studios, Zurich who made this course possible. It was a great learning experience and got me really excited about computer graphics.
The website template was borrowed from Michaël Gharbi.

Computer Graphics 2022 ETH Zurich Matias Turkulainen

Introduction

Motivation

Simple Features

Lens distortion

Depth of Field

Image Textures

Procedural Textures

Rendering on a High Performance Cluster

Advanced Features

Environmental Map Emitter

Implementation Details

Evaluation

Homogeneous Participating Media

Implementation Details

Evaluation

Final Render

Acknowledgements

Computer Graphics 2022
ETH Zurich
Matias Turkulainen