Iceman Project

This is the final rendering competition for CS87/287 Fall Rendering Algorithms. The course focuses on physically based rendering using Monte Carlo methods. This year's theme is "Reflections on Time." During the last two weeks of the term, students individually implement advanced rendering algorithms based on Darts, a C++17 educational Monte Carlo ray-tracer skeleton, and produce a photorealistic image with their renderer.

For inspiration, I drew from Nele Azevedo's installation "Melting Men."^[1]Her work elegantly evokes the passage of time: the figures gradually melt and disappear. Borrowing that motif, I created an "Iceman" whose head droops slightly, suggesting contemplation. Facing him is a fire, providing a strong visual and semantic contrast.

On the technical side, I implemented volume path tracing using the null-scattering framework^[2][3]. I also added importance sampling for the environment map.

Heterogeneous Media for Realistic Ice

To make the ice look realistic, I observed that natural ice is not perfectly transparent. The interior often appears foggy. As water freezes, dissolved air is expelled and becomes trapped as bubbles, producing a cloudy center. To reproduce this, I render the interior as a volume whose density increases toward the center of the geometry. This required three subtasks: define the volume boundary, set the volume density, and render the volume. I handled the first two tasks in Blender using its fluid simulation. By setting emission sources, boundaries, and collisions, the simulated fluid gets trapped inside the Iceman mesh. The cached fluid data was located at:
C:\Users\MyName\AppData\Local\Temp\blender_a34288\cache_fluid_4b85f272\data
By tweaking the simulation parameters we can create an interior volume inside the ice. (Note: fluid tasks are often easier and more powerful in Houdini though.)

I modeled the figure in Blender and obtained background assets such as trees, an environment texture, and fire.vdb from the internet.

To render the heterogeneous volume, I used the null-scattering framework. The idea is simple and straightforward to implement. Recall the volume rendering equation (assuming no surfaces for simplicity):

• $L_i(\mathbf{x},\vec{\omega})$ — Incoming radiance at point x in direction ω.
• $f_t(\mathbf{x},\mathbf{x}_t)$ — Probability density of light colliding with a particle at the location. Also called free-flight PDF. This PDF is only determined by the total "density" of the particles, $\sigma_t(\mathbf{x}_t)$, at location $\mathbf{x}_t$.
• $\frac{\sigma_a(\mathbf{x}_t)}{\sigma_t(\mathbf{x}_t)}$ — Proportion of emitting particles at the location. / Probability of interacts with an emitting particle at the location.
  There's a common confusion here: the $\sigma_a$should be the density of absorption particles, not emitting one. We can introduce another density for absorption for physical correctness, however, in most cases, for simplicity we use absorption coefficient times Le term. This means the Le term here is not physically correct. But it is totally fine for artistic control.
• $\frac{\sigma_s(\mathbf{x}_t)}{\sigma_t(\mathbf{x}_t)}$ — Proportion of scattering particles at the location. / Probability of interacts with a scattering particle at the location.
• $L_e$ — Emitted radiance from sources inside the medium. The clarification above fits here.
• $L_s$ — Scattered radiance (in-scattered light).

In a homogeneous medium the total extinction $\sigma_t(\mathbf{x}_t)$ is constant, making sampling trivial (sampling a 1D free-flight PDF). Null-scattering introduces "null particles" to homogenize a heterogeneous medium: when a null particle is sampled, the ray continues unchanged. The transport equation becomes:

You then sample travel distances according to the free-flight PDF defined by the (homogenized) medium. Intuitively, null-scattering reduces the complexity of sampling in heterogeneous media by restricting the free-flight PDF to a 1D function. This enables rendering the fog efficiently.

Importance Sampling the Environment Map

The environment map is a spherical light source at infinite radius represented by a 2D image. There are several mappings between the sphere and the plane, e.g. spherical coordinates, octahedral mapping, and equal-area mapping. I used spherical coordinates, whose drawback is non-uniform area sampling: For each altitude angle $\theta$, there are equal numbers of $\phi$ texels in the 2D image. However, the area of each band within infinitesimal altitude $\theta \quad \textit{to} \quad \theta+ \Delta\theta$ gets smaller as $\theta$ increases. Therefore, the precision of the environment map at the pole will be higher than the equator.

The conversion between solid angle and spherical coordinates is:$dA = \sin\theta\,d\theta\,d\phi$

The light transport equation in spherical coordinates is:

We importance-sample the sine-weighted radiance to reduce variance. Practically, we precompute a sine-weighted texture and sample it as a 2D piecewise-constant PDF. The sine-weighted texture looks like this:

A detail to notice when implementing this: since we are sampling a 2D piecewise constant PDF within a rectangle, there's a 2-step mapping to get the correct PDF for spherical coordinate: from rectangle to a unit square, then to spherical coordinate. Need to consider 2 Jacobians for change of measure. Related code is down here:

float EnvMaterial::pdf(const Ray3f& ray, const Vec3f& scattered, const HitRecord& hit) const {
                    if (m_distribution2D == nullptr) {
                        return INV_FOURPI; // uniform over sphere/hemisphere depending on use
                    }
                    Vec3f dir = normalize(scattered);
                    // Convert to env-map convention if needed (flip z)
                    Vec3f dir_rh_tocenter(dir.x, dir.y, -dir.z);
                    // Map direction to equirectangular UV where u,v in [0,1)
                    Vec2f uv = Spherical::direction_to_equirectangular(dir_rh_tocenter);
                    // Change-of-variables from (u,v) to solid angle ω: dω = sin(θ) dθ dφ with θ = vπ, φ = u2π
                    // |J| = sin(θ) * 2π^2 => pdf_ω = pdf_uv / (sin(θ) * 2π^2)
                    float pdf = m_distribution2D->pdf(uv) / (2.0f * M_PI * M_PI * sin(uv.y * M_PI));
                    return pdf;
                }