Nori ↑

Introduction ↑

Building upon the Nori renderer from ETH Zurich's Computer Graphics course (FS2025), Nori represents my exploration into advanced rendering techniques. Here I document the features I implemented beyond the course requirements, along with my notes and validations.

Emitters ↑

Spot Light ↑

Spotlights emit light in a cone of directions from their position. This is discussed in chapter 12.2.1 of the PBR book.

Description ↑

Files created or modified:

• src/spotlight.cpp

Validation scenes:

• scenes/nori-beyond/spotlight/cbox.xml
• scenes/nori-beyond/spotlight/cbox_mitsuba.xml

A new emitter type SpotLight is implemented in the file spotlight.cpp.

The origin of the light is at $(0, 0, 0)$, and points towards the $+z$ axis. A spot light is defined by two angles -- falloffStart, where everything inside this angle is fully illuminated (intensity multiplier = 1), and totalWidth, where illumination fades to black at this angle, and everything outside is dark (intensity multiplier = 0). To avoid expensive $\arccos$ operations, I store the the cosines of the angles cosFalloffStart and cosTotalWidth, and since $\cos \theta$ decreases as $\theta$ increases, $\cos \theta_\text{fall off start} > \cos \theta_\text{total width}$. To create a soft transition zone (penumbra) in the outer cone, I use the SmoothStep function, which is a cubic function (Hermite interpolation), that takes in $\cos \theta_\omega, \cos \theta_\text{total width}, \cos \theta_\text{fall off start}$, where $\cos \theta_\omega$ is the cosine angle between $+z$ and the vector with direction $\mathbf{\omega}$ pointing from the light origin to the point in scene:

1. Compute $t = \frac{\cos \theta_\omega - \cos \theta_\text{total width}}{\cos \theta_\text{fall off start} - \cos \theta_\text{total width}}$ and $t = 0$ if $\cos \theta_\omega < \cos \theta_\text{total width}$ or $\cos \theta_\omega > \cos \theta_\text{fall off start}$.
2. Return $t^2(3 - 2t)$ (Hermite interpolation)

The angular intensity $I(\mathbf{\omega})$ is the radiant intensity in the direction $\mathbf{\omega}$ (power emitted per solid angle in direction $\mathbf{\omega}$), and is given by

$$I(\mathbf{\omega}) = \underbrace{\text{scale} \times \text{spectrum}}_\text{Peak intensity in the forward direction before angular attenuation} \times \underbrace{\text{SmoothStep}(\cos \theta_\omega, \cos \theta_\text{total width}, \cos \theta_\text{fall off start})}_\text{Scalar in $[0, 1]$ depending on $\theta_\omega$}$$

A surface patch $\mathrm{d}A$ at point $\mathbf{p}$ subtends a small solid angle $\mathrm{d}\mathbf{\omega}$ around direction $\mathbf{\omega}$, with the relation $\mathrm{d}\mathbf{\omega} = \frac{\cos \theta}{r^2}\mathrm{d}A$, where $\theta$ is the angle between the surface normal at $\mathbf{p}$ and direction from $\mathbf{p}$ to the light. The power arriving at $\mathrm{d}A$ from that infinitesimal solid angle is given by

$$\mathrm{d}\Phi = I(\mathbf{\omega})\mathrm{d}\mathbf{\omega} = I(\mathbf{\omega})\frac{\cos \theta}{r^2}\mathrm{d}A$$

The irradiance on the patch is $E = \frac{\mathrm{d}\Phi}{\mathrm{d}A} = \frac{I(\mathbf{\omega})\cos \theta}{r^2}$. The radiance $L$ is defined so that

$$\begin{align*} \mathrm{d}\Phi &= L \cos \theta \mathrm{d}A\mathrm{d}\mathbf{\omega} \\ L &= \frac{\mathrm{d}\Phi}{\cos \theta \mathrm{d}A\mathrm{d}\mathbf{\omega}} \\ &= \frac{I(\mathbf{\omega})\mathrm{d}\mathbf{\omega}}{\cos \theta \mathrm{d}A\mathrm{d}\mathbf{\omega}} \\ &= \frac{I(\mathbf{\omega})}{r^2} \end{align*}$$

Derivation of flux (not used explicitly in this implementation because of the uniform sampling strategy)

A spherical cap (cone of half-angle $\theta$) has solid angle $\Omega(\theta) = 2\pi(1 - \cos \theta)$. Within the region of falloffStart, intensity is constant, so flux in that region is

$$\Phi_\text{fall off start} = I_\text{fall off start}\Omega(\theta_\text{fall off start}) = 2\pi I_\text{fall off start}(1 - \cos \theta_\text{fall off start})$$

Between falloffStart and totalWidth, the intensity decreases smoothly to 0. The PBR book approximates the integral by assuming that intensity in the region averages to half the peak value. The solid angle of the penumbra annulus is $\Omega_\text{penumbra} = 2\pi(\cos \theta_\text{fall off start} - \cos \theta_\text{total width})$. If the average intensity is $I_\text{fall off start} / 2$, then

$$\Phi_\text{penumbra} \approx \frac{I_\text{fall off start}}{2}\Omega_\text{penumbra} = = 2\pi I_\text{fall off start}\frac{\cos \theta_\text{fall off start} - \cos \theta_\text{total width}}{2}$$

Total flux is thus

$$\Phi \approx \Phi_\text{fall off start} + \Phi_\text{penumbra} = 2\pi I_\text{fall off start}\left[(1 - \cos \theta_\text{fall off start}) + \frac{\cos \theta_\text{fall off start} - \cos \theta_\text{total width}}{2}\right]$$

SpotLight::eval(...)

Delta emitters have zero area so there is measure-zero probability of hitting them via BSDF sampling.

SpotLight::sample(...)

First assign the fields of the emitter query record lRec given the position of the light and the reference point $\mathbf{p}$. Then compute the angle attenuation term using $\text{SmoothStep}$ given the cosine of the angle between the m_direction (direction of spot light) and the vector pointing from the light to $\mathbf{p}$ (lRec.ref). Finally compute the intensity and return the radiance.

SpotLight::pdf(...)

Delta lights have zero measure, thus PDF = 0. This is required for MIS.

Spotlight::samplePhoton(...)

First sample a direction from the light using squareToUniformSphereCap and assign the sampled direction to ray. Then compute the angle attenuation term using $\text{SmoothStep}$ given the cosine of the angle between the m_direction (direction of spot light) and the sampled direction. Finally return the photon power, which is the intensity divided by the sampling PDF.

Validation ↑

I compared my implementation with Mitsuba with the Cornell box scene with two cubes, and a spot light positioned behind the right one, pointing towards the left one.

I also updated path_mis.cpp to support delta lights (spotlights and point lights) by modifying the MIS weight calculation. For delta lights (detected by checking if ems_lRec.n.isZero()), I set the MIS weight to 1.0 since delta lights cannot be sampled by BSDF sampling (the probability of randomly hitting a point light via BSDF is zero). This ensures that only emitter sampling contributes when rendering scenes with delta lights, preventing the image from being incorrectly darkened by MIS weighting.

There is a reasonable difference in rendered images. In my implementation, after computing $t$, I apply a Hermite interpolation (SmoothStep) following the PBR book, which creates an "S-curve" falloff. My spotlight hence stays brighter for longer near the center (the "shoulder" of the S-curve) and fade out more gradually near the absolute edge (the "toe" of the S-curve) compared to Mitsuba's sharp linear drop, which applies applies a pure linear ram. Also, my interpolation domain is cosine where Mitsuba is the angle itself. I believe that the visual differences make sense based on these two implementation differences.

Directional Light ↑

Directional Light is also known as distant light, which desccribes an emitter that deposits illumination from the same direction at every point in space. This feature is discussed in chapter 12.3 of the PBR book.

Description ↑

Files created or modified:

• src/directionallight.cpp

Validation scenes:

• scenes/nori-beyond/directional_light/

A new emitter type DirectionalLight is implemented in the file directionallight.cpp.

The radiance at a point $\mathbf{x}$ from direction $\boldsymbol{\omega}_\text{i}$ due to directional light with direction $\boldsymbol{\omega}_\ell$ is given by

$$L_\text{i}(\mathbf{x}, \boldsymbol{\omega}_\text{i}) = L_\text{e}\delta(\boldsymbol{\omega}_\text{i} - \boldsymbol{\omega}_\ell)$$

The physical reason that directional lights make sense is that for point lights, $L_\text{i} \propto \frac{1}{r^2}$ where $r$ is the distance from the point light. For directional lights, $r \to \infty$. However, the emitting area also grows as $r^2$. The two effects cancel.

DirectionalLight::eval(...)

Delta emitters have zero area so there is measure-zero probability of hitting them via BSDF sampling.

DirectionalLight::sample(...)

Given a shading point $\mathbf{x}$, we sample the incoming direction $\boldsymbol{\omega}_\text{i}$ deterministically since there is only one valid direction. We set:

• $\boldsymbol{\omega}_\text{i} = -\boldsymbol{\omega}_\ell$ (toward the light source)
• $\mathbf{p} = \mathbf{x} + \boldsymbol{\omega}_\text{i} \cdot 2r_\text{scene}$ (a fictitious point far away)
• $\text{pdf} = 1$ (delta distribution, handled specially)
• Shadow ray from $\mathbf{x}$ toward $\boldsymbol{\omega}_\text{i}$ with infinite extent

The returned radiance is simply $L_\text{e} = \texttt{scale} \times \texttt{spectrum}$. Unlike area lights, there is no $1/r^2$ falloff since the light is infinitely distant.

DirectionalLight::pdf(...)

Delta lights have zero measure, thus PDF = 0. This is required for MIS to correctly handle the delta distribution---the MIS weight computation uses this PDF, and returning 0 ensures that BSDF-sampled paths that happen to align with the light direction don't incorrectly contribute.

DirectionalLight::samplePhoton(...)

For photon mapping, we need to emit photons that cover the entire scene uniformly from the light direction. Following PBR book 12.5, we:

1. Sample a point on a disk of radius $r_\text{scene}$ perpendicular to the light direction $\boldsymbol{\omega}_\ell$
2. Position this disk behind the scene (offset by $r_\text{scene}$ along $-\boldsymbol{\omega}_\ell$)
3. Emit a ray from this point in direction $\boldsymbol{\omega}_\ell$

The disk sampling uses an orthonormal frame built from $\boldsymbol{\omega}_\ell$:

$$\mathbf{o} = r_\text{scene}(u \cdot \mathbf{s} + v \cdot \mathbf{t}) - r_\text{scene} \cdot \boldsymbol{\omega}_\ell$$

where $(u, v)$ is a uniform disk sample and $(\mathbf{s}, \mathbf{t}, \boldsymbol{\omega}_\ell)$ form an orthonormal basis.

The photon power is $L_\text{e} \times \pi r_\text{scene}^2$ (radiance times the disk area), which accounts for the uniform PDF of $1/(\pi r_\text{scene}^2)$ over the sampling disk.

Validation ↑

Environment Map Emitter ↑

An environment map emitter (infinite area light) wraps the scene in an infinitely large sphere textured with an HDR image, so that every direction in space receives radiance defined by the map. This is discussed in chapter 12.4 of the PBR book.

Description ↑

Files created or modified:

• src/envmap.cpp
• include/nori/emitter.h
• include/nori/scene.h and src/scene.cpp
• All integrators: src/path_mis.cpp, src/path_mats.cpp, src/path_ris.cpp, src/direct_mis.cpp, src/direct_mats.cpp, src/direct_ems.cpp, src/direct.cpp
• envmap_sampling_viz.py
• latlong_to_equalarea.py

Validation scenes:

• scenes/nori-beyond/environment_map/cbox_envmap.xml

A new emitter type EnvironmentMap is implemented in the file envmap.cpp, registered as "envmap".

For every direction $\boldsymbol{\omega}$ in the sphere, light comes from that direction with some radiance $L(\boldsymbol{\omega})$. If we shoot a ray and it misses everything in the scene, we look up the ray direction in the HDR environment image. However, since our HDR file is a 2D image, but directions live on a sphere, we need a mapping

$$\boldsymbol{\omega} = (x, y, z) \leftrightarrow (u, v) \in [0, 1]^2$$

Following Chapter 3.8.3 in the PBR book, the environment map image uses the equal-area octahedral parameterization rather than the traditional equirectangular (latitude-longitude) projection. The traditional equirectangular mapping $u = \phi / 2\pi, v = \theta / \pi$ compresses solid angles near the poles ($\sin\theta \to 0$), so polar pixels cover vanishingly small solid angles. This requires an explicit $\sin\theta$ correction in the sampling distribution and introduces a $1/(2\pi^2 \sin\theta)$ Jacobian that becomes singular at the poles.

The equal-area octahedral mapping avoids these issues. It combines:

1. Octahedral projection: Project the unit sphere onto the faces of an octahedron, then unfold the octahedron into a square
2. Concentric disk-to-hemisphere mapping: Apply an area-preserving mapping (from Shirley & Chiu 1997) so that each pixel subtends exactly the same solid angle: $4\pi / (W \times H)$

Forward mapping (sphere to square, implemented in src/envmap.cpp as a helper function equalAreaSphereToSquare(const Vector3f &d)): Given a direction $\boldsymbol{\omega} = (x, y, z)$ on the unit sphere, the mapping computes $(u, v) \in [0, 1]^2$. The setup is as follows: Imagine an octahedron inscribed in a sphere. The ocatahedron can be unfolded to a square -- the top four faces of it (where $z \geq 0$) map directly to a diamond in the center of the 2D square, and the bottom four faces (where $z < 0$) are "folded out" into the four corner triangles.

1. In the setup described above, we can then define the $\mathrm{L}1$ radius $R$ from latitude. The latitude (circle on the sphere at a specific height $z$) is mapped to the square $[-1, 1]^2$ such that the point on the square sits on a concentric diamond around the origin if $z \geq 0$ or a concentric corner triangle near $(\pm 1, \pm1)$ if $z < 0$. This is defined by $R = \sqrt{1 - |z|}$, which means that:
   • when $z \geq 0$, the $\mathrm{L}1$ radius on the square is given by $|u| + |v| = R$
   • when $z < 0$, the $\mathrm{L}1$ radius on the square is given by $|u| + |v| = 2 - R$
2. Next we determine where along that line the point lands. We first work in the positive quadrant of $[-1, 1]^2$ using $|x|$ and $|y|$. Let

   double a = max(fabs(x), fabs(y));
   double b = min(fabs(x), fabs(y));
   b = (a == 0) ? 0 : b / a;
   

   Then we approximate $\phi = \frac{4}{\pi}\arctan(b)$ using a minimax polynomial. In implementation, the exact coefficients in Nori are taken from PBRT's source code, in src/pbrt/util/math.cpp inside EqualAreaSphereToSquare(), where the authors mention that the code is "Via source code from Clarberg: Fast Equal-Area Mapping of the (Hemi)Sphere using SIMD". This avoids calling atan() directly:

   double phi = P(b);   // approximates (4/pi) * atan(b)
   if (fabs(x) < fabs(y))
       phi = 1 - phi;
   

   This gives a normalized position along the diamond edge in the positive quadrant. The corresponding coordinates are

   double vp = phi * R;
   double up = R - vp;
   

   For the upper hemisphere ($z \ge 0$), this already lies on the inner diamond $|u| + |v| = R$. For the lower hemisphere ($z < 0$), the point must be folded out to the corner triangle:

   if (z < 0) {
       swap(up, vp);
       up = 1 - up;
       vp = 1 - vp;
   }
   

   After this step, $(u_p, v_p)$ lies in the positive quadrant with the correct $\mathrm{L}1$ radius: $R$ for the upper hemisphere and $2 - R$ for the lower hemisphere.
3. We then restore the original $x$ and $y$ signs and remap from $[-1, 1]^2$ to $[0, 1]^2$:

   double u_signed = copysign(up, x);
   double v_signed = copysign(vp, y);
   
   double u = (u_signed + 1) / 2;
   double v = (v_signed + 1) / 2;
   

   This places the point in the correct quadrant of the final square.

Inverse mapping (square to sphere, implemented in src/envmap.cpp as a helper function equalAreaSquareToSphere(const Point2f &p)): Given $(u, v) \in [0, 1]^2$, we recover $\boldsymbol{\omega}$. Following PBRT, it uses a direct branch-free formula that maps the square back to the sphere.

1. Remap $(u, v)$ from $[0, 1]^2$ to $[-1, 1]^2$ and work with $(u', v')$ in that signed square. Their absolute values determine how far the point is from the diagonal $|u'| + |v'| = 1$, which separates the upper and lower hemispheres:

   double up = fabs(u');
   double vp = fabs(v');
   double signedDistance = 1 - (up + vp);
   double r = 1 - fabs(signedDistance);
   

   Here, signedDistance determines the hemisphere: positive means upper hemisphere, negative means lower hemisphere. The quantity $r$ is the same equal-area radius used by the forward map.

2. Reconstruct the angular parameter within the positive quadrant:

   double phi = (r == 0 ? 1 : (vp - up) / r + 1) * (M_PI / 4);
   

   This is the inverse of the forward-map relation that used $\phi$ to split the radius into two square coordinates. It recovers the angle within one octahedral sector.

3. Convert $(r, \phi)$ back to a 3D direction:

   double z = copysign(1 - r * r, signedDistance);
   double cosPhi = copysign(cos(phi), u');
   double sinPhi = copysign(sin(phi), v');
   double rFactor = r * sqrt(2 - r * r);
   return Vector3f(cosPhi * rFactor, sinPhi * rFactor, z);
   

   These steps are not direct inverses of the swap() / 1 - u folding logic used in equalAreaSphereToSquare(). Instead, they reconstruct the sphere point from the equal-area square parameterization in one shot:

   1. z = copysign(1 - r^2, signedDistance) recovers the hemisphere and height on the sphere.
   2. cosPhi and sinPhi recover the azimuthal direction inside the correct quadrant by attaching the signs of $u'$ and $v'$.
   3. rFactor = r\sqrt{2-r^2} is the radius of the projected point in the $xy$-plane, coming from PBRT's equal-area disk-to-hemisphere mapping.

To find the light coming from direction $\boldsymbol{\omega}$, we convert that direction into map coordinates $(u, v)$, find which pixel $(x, y)$ that lands on in our $W \times H$ image, take that pixel's color and multiply it by a brightness scale:

$$L(\boldsymbol{\omega}) = \texttt{scale} \times \texttt{image}\left[\lfloor v \cdot H \rfloor\right]\left[\lfloor u \cdot W \rfloor\right]$$

Uniform sphere sampling would be very inefficient since most of the radiance is typically concentrated in a small region (e.g. the sun). Following PBRT Chapter 12.4, we build a 2D piecewise-constant distribution over the image pixels.

Because the equal-area mapping ensures every pixel covers the same solid angle, the sampling weights are simply the luminance of each pixel, and no $\sin\theta$ correction is needed. This is one of the key advantages of the equal-area mapping.

The distribution is built in two levels:

1. Conditional CDFs (per row): For each row $y$, build a CDF over columns weighted by $w(x, y) = \text{luminance}(x, y)$
2. Marginal CDF (over rows): Build a CDF over rows weighted by the integral of each row $\sum_x w(x, y)$

Sampling draws a row from the marginal CDF, then a column from that row's conditional CDF, both via binary search. The sampled pixel $(x, y)$ is converted back to a direction via EqualAreaSquareToSphere.

The PDF in solid-angle measure is obtained by a simple change of variables. Since the equal-area mapping has Jacobian $1/(4\pi)$ (mapping $[0, 1]^2$ to the unit sphere of area $4\pi$):

$$p(\boldsymbol{\omega}) = \frac{p(u, v)}{4\pi}$$

where $p(u, v) = p_\text{marginal}(v) \times p_\text{conditional}(u \mid v) \times W \times H$ is the discrete probability density in image space.

include/nori/emitter.h

Two functions isEnvironment and evalEnvironment are added. isEnvironment returns true if this is an environment emitter and is false by default. EnvironmentMap::isEnvironment(...) overrides to true. evalEnvironment is needed for the case where a ray escapes the scene entirely. With no intersection information, we only have the ray's direction. evalEnvironment is a direction-only convenience method added to this base Emitter class for this specific case.

Scene::addChild(...)

In the header file, we add a single stored pointer to the environment light m_envEmitter and the method getEnvironmentEmitter that returns the pointer to the environment emitter.

In the case of an emitter, we additionally check if the emitter is an environment light and if pointer m_envEmitter is already non-null, i.e. the scene already has an environment emitter stored, in which case we throw an exception.

EnvironmentMap::eval(...)

Implements EnvironmentMap::evalEnvironment(...) with the direction lRec.wi. The evaluation proceeds in the following steps:

1. Transform the world-space direction $\mathbf{d}$ on the unit sphere into 2D texture coordinates $(u, v) \in [0, 1]^2$ using equalAreaSphereToSquare

2. Convert the continuous $(u, v)$ coordinates to pixel space with pixel-centered addressing:

   $$f_x = u \cdot W - 0.5, \quad f_y = v \cdot H - 0.5$$

   The $-0.5$ offset treats pixel centers as integer coordinates.

3. Compute the integer coordinates of the four corners:

   $$x_0 = \lfloor f_x \rfloor, \quad y_0 = \lfloor f_y \rfloor, \quad x_1 = x_0 + 1, \quad y_1 = y_0 + 1$$

   Compute the fractional offsets:

   $$t_x = f_x - x_0, \quad t_y = f_y - y_0$$

   Clamp all coordinates to the valid range $[0, W-1]$ and $[0, H-1]$ respectively to handle boundary cases.

4. Fetch the four corner pixel values from m_data:

   $$c_{00} = \texttt{m\_data}[y_0 \cdot W + x_0], \quad c_{10} = \texttt{m\_data}[y_0 \cdot W + x_1]$$ $$c_{01} = \texttt{m\_data}[y_1 \cdot W + x_0], \quad c_{11} = \texttt{m\_data}[y_1 \cdot W + x_1]$$

   Blend using bilinear interpolation:

   $$L(\mathbf{d}) = \texttt{scale} \times \left[(1 - t_x)(1 - t_y) \, c_{00} + t_x(1 - t_y) \, c_{10} + (1 - t_x)t_y \, c_{01} + t_x t_y \, c_{11}\right]$$

   This provides smooth filtering across pixel boundaries, reducing aliasing compared to nearest-neighbor lookup.

buildDistribution()

This is a method that fills in EnvironmentMap's member variables m_marginalPdf, m_marginalCdf, m_conditionalPdf and m_conditionalCdf using m_width, m_height, m_data. m_conditionalCdf is an array of size m_height * (m_width + 1) since each row requires (m_width + 1) boundary points to define its m_width bins and there are m_height rows in total. m_marginalCdf is an array of size m_height + 1 for row sampling, according to a marginalWeights array.

It follows the steps:

1. Build the conditional CDF for each row. For each row $y$, compute the base (start of that row's CDF inside the flat m_conditionalPdf array) and initialize m_conditionalPdf[base] = 0.
2. For each pixel $(x, y)$ in row $y$, retrieve the luminance by looking up the pixel value in m_data. Then build a prefix sum asnd store that into m_conditionalCdf[base + x + 1].
3. The total brightness of the row is then stored in m_conditionalCdf[base + m_width], which is then saved in marginalWeights[y] to inform row sampling (brighter rows should be sampled more often).
4. If the row has non-zero brightness, normalize the row CDF so that it becomes a proper CDF in $[0, 1]$. We also fill in m_conditionalPdf using the luminance of each pixel in the row normalized by this row total. Otherwise, it means that the row is completely black, and to avoid division by zero, we fall back to uniform sampling, where each column of the row gets probability 1 / m_width.
5. After building m_conditionalCdf and m_conditionalPdf, we build m_marginalCdf and m_marginalPdf using marginalWeights. Similar to the previous step, we normalize m_marginalCdf and m_marginalPdf.

buildDistribution is called once after the image data is loaded. The CDFs m_marginalCdf and m_conditionalCdf are then used to sample a pixel on the environment map, and m_marginalPdf and m_conditionalPdf are used to compute the Monte Carlo radiance of the sampled pixel.

EnvironmentMap::sample(...)

1. Sample a row $v$ from the marginal CDF using the first random number and record $p_\text{marginal}(v)$
2. Sample a column from that row's conditional CDF using the second random number and and record $p_\text{marginal}(u | v)$
3. Convert the sampled pixel to the coordinates in a $[0, 1]^2$, then to a direction $\boldsymbol{\omega}$ using EqualAreaSquareToSphere
4. Compute the solid-angle PDF $p(\boldsymbol{\omega}) = p(u, v) / (4\pi)$ where $p(u, v) = p_\text{marginal}(v) p_\text{marginal}(u | v)$
5. Set lRec.p to a distant point along $\boldsymbol{\omega}$ (at $2 r_\text{scene}$), and the shadow ray to $[\mathbf{x}, \boldsymbol{\omega}, \epsilon, \infty)$
6. Return $L(\boldsymbol{\omega}) / p(\boldsymbol{\omega})$.

EnvironmentMap::pdf(...)

We convert the direction lRec.wi to image coordinates $(u, v)$ using EqualAreaSphereToSquare and then look up the marginal and conditional PDFs, and returns $p(u, v) / (4\pi)$.

EnvironmentMap::samplePhoton(...)

For photon mapping, we sample a direction from the environment map's importance distribution, then place the photon origin on a disk of radius $r_\text{scene}$ perpendicular to that direction (similar to the directional light approach). The photon power is $L(\boldsymbol{\omega}) \times \pi r_\text{scene}^2$.

Integrator Modifications

The environment map is a non-geometric emitter so rayIntersect can never find it. Two changes were needed in every integrator:

1. Camera rays that miss geometry: Instead of returning black, check for an environment emitter and return $L(\boldsymbol{\omega}_\text{ray})$. This renders the environment as the background.
2. BSDF-sampled rays that escape the scene: Evaluate the environment map and add its contribution. For MIS integrators, the MIS weight is computed using the environment map's solid-angle PDF as the emitter PDF.

Emitter sampling already works automatically since the environment map is added to the scene's emitter list and implements sample() and pdf().

Validation ↑

The environment used for validation here was taken from PolyHaven:

which was converted to an equal area map using latlong_to_equalarea.py. The upper hemisphere (the sky) is mapped to the centre diamond, and the lower hemisphere (the ground) is mapped to the four corners, just as described in the description above. Sampling is proportional to the luminance, as seen in the heat map below (brighter pixels represent more frequent sampling).

The validation scene against Mitsuba:

Textures ↑

Images as Texture ↑

Description ↑

Files created or modified:

• include/nori/lodepng.h: downloaded from lodepng
• src/lodepng.cpp: downloaded from lodepng
• src/imagetexture.cpp

Validation scenes:

• scenes/nori-beyond/image_texture/mesh-image-texture.xml
• scenes/nori-beyond/image_texture/mesh-image-texture-mitsuba.xml

This feature enables the use of PNG images as textures to replace constant albedo colors on objects in the scene. A new texture type ImageTexture is implemented in imagetexture.cpp. The texture supports tiling and scaling.

The standard texture mapping pipeline includes scene definition -> texture loading at startup -> ray intersection -> UV coordinate interpolation -> texture lookup -> colour / material evaluation. To breakdown the implementation, the following details the pipeline:

1. Image texture loading at startup: An image loader would be needed to parse the image file in order to handle decompressing data, convert to RGB format, etc.
2. Ray intersection: When a ray hits a mesh triangle, I would need to get the barycentric coordinates of the hit point within the triangle bary_coordinates.
3. UV coordinate interpolation: The 3D mesh has UV coordinates stored per vertex, normalized to $[0, 1]$. For example, Point2f uv0 = m_UV.col(v0) is the UV coordinates at vertex 0. Using bary_coordinates, I can interpolate the UV coordinates: its.uv = bary_coordinates.x() * uv0 + bary_coordinates.y() * uv1 + bary_coordinates.y() * uv2. This part and the previous two parts are already handled in nori.
4. Texture lookup: I need to loopk up the colour in the loaded texture data. This consists of first converting the UV coordinates ($[0, 1]$) to pixel coordinates ([width, height]), then looking up the image loader data.
5. Colour/Material evaluation: This is implemented in the BSDF eval methods (e.g. Diffuse's eval method calls m_albedo->eval(bRec.uv)). Finally, the BSDF value is evaluated with bsdf->eval(bRec).

For the image loader, I searched for commonly used ones and were deciding between lodepng and stb_image. I decided to use lodepng since I only need PNG support and that it has separate compilation and has automatic memory management.

For the PNG image loading, I first construct the flattened 1D array m_data in which each entry is a Color3f, for memory efficiency.

ImageTexture::eval(...)

Given the UV coordinates on the 3D surface, this method should find the corresponding pixel in the image texture, and return its RGB colour.

First, apply the tiling and scaling by computing the fractional part of the scaled UV coordinates

float img_u = std::fmod(uv.x() * m_scale.x(), 1.0f);
float img_v = std::fmod(uv.y() * m_scale.y(), 1.0f);

For example, if uv = (0.75, 0.75) and m_scale = (2.0, 2.0), then after scaling, I get (1.5, 1.5) which is beyond $[0, 1]$ and I are in the second repetition, so I subtract by $\lfloor 1.5 \rfloor = 1$ to get $(0.5, 0.5)$, which means that at surface point $(0.75, 0.75)$ I sample texture from $(0.5, 0.5)$, because the texture is repeated twice.

Next, convert the UV coordinates to pixel coordinates by scaling img_u by m_width and img_v by m_height and trancating to integer pixel index. For the V coordinate, I need to flip the axis because in the image storage, row 0 is at the top of the image but V = 0 is at the bottom of the texture. Then I clamp the bounds to prevent the pixel coordinates being negative or going beyond the last column or row.

Finally, convert from 2D to 1D array index to look up the RGB value from m_data which I constructed.

Validation ↑

For the validation scenes, I reuse the scene mesh-texture from assignment 1

Normal mapping ↑

Normal mapping replaces the shading normal with a normal sampled from a texture (RGB image). The texture stores perturbations of the normal direction. This feature is discussed in chapter 10.5.3 of the PBR book.

Description ↑

Files created or modified:

• src/normalmap.cpp
• include/nori/bsdf.h
• src/diffuse.cpp
• src/mesh.cpp

Validation scenes:

• scenes/nori-beyond/normal_mapping/mesh-normal-map.xml
• scenes/nori-beyond/normal_mapping/mesh-no-normal-map.xml
• scenes/nori-beyond/normal_mapping/mesh-normal-map-mitsuba.xml
• scenes/nori-beyond/normal_mapping/mesh-no-normal-map-mitsuba.xml

The current rendering pipeline follows ray generation -> ray-scene interaction -> mesh.cpp::setHitInformation which fills in its -> integrators to compute radiance. The NormalMap texture would return the Normal3f at a given UV coordinate in the 3D surface. We need to modify mesh.cpp::setHitInformation so that its.shFrame is the perturbed frame according to the normal map.

src/normalmap.cpp

A new texture NormalMap is implemented in normalmap.cpp. Normal maps are stored in the tangent space (aligned with per-vertex tangent, bitangent, normal) using (R, G, B) = (nx, ny, nz) with each component in $[0, 1]$. Since real normals are in $[-1, 1]$, the PBR book decodes them by x = 2 * r - 1, and correspondingly for y and z. After that, we store the Normal3f(x, y, z) into m_data. In include/nori/frame.h, the Frame class defines s and t as tangent vectors and n as the normal and the key transformation method is:

Vector3f toWorld(const Vector3f &v) const {
    return s * v.x() + t * v.y() + n * v.z();
}

which is a standard right-handed coordinate frame with no special tangent convention defined. So for normal maps with OpenGL format, we pass through (x, y, z) which will be mapped to (s, t, n) and for normal maps with DirectX format, we pass (x, -y, z) since DirectX conventions have Y pointing down. Other than this detail and that the final returned Normal3f should be normalized, the implementation of eval is the same as as that in ImageTexture.

include/nori/bsdf.h

We added the method getNormalMap() which returns true if the BSDF has a normal map. We need to add this method because image textures are stored in m_albedo within the BSDF (diffuse.cpp::addChild()), and then the sample method uses it. However, since we need to modify its.shFrame inside mesh.cpp::setHitInformation, it must access m_normalMap.

src/diffuse.cpp

To use the normal map, we need at least one BSDF that stores it. We use the diffuse BSDF. For each instance where m_albedo appears, we also add m_normalMap. This includes the constructor, destructor, addChild() and toString(). In addition, we also implement the getNormalMap() override to return the normal map so the mesh can retrieve it.

Mesh::setHitInformation()

After interpolating base normal from mesh, we need to query BSDF for the normal map using getNormalMap(). If normal map exists, we evaluate the normal map at its.uv to get the tangent-space normal. Then we compute the tangent frame from UV gradients and transform the normal to world space. Then we update its.shFrame with the perturbed normal. In detail, we first read the normal map $\mathbf{n}_s = (n_x, n_y, n_z)$ using Normal3f mappedNormal = normalMap->eval(its.uv), where eval was implemented in normalmap.cpp.

Let the triangle vertex positions be $\mathbf{p}_0, \mathbf{p}_1, \mathbf{p}_2 \in \mathbb{R}^3$ and the triangle UVs be $\mathbf{u}_0, \mathbf{u}_1, \mathbf{u}_2 \in \mathbb{R}^2$. We define $\Delta \mathbf{p}_1 = \mathbf{p}_1 - \mathbf{p}_0$, $\Delta \mathbf{p}_2 = \mathbf{p}_2 - \mathbf{p}_0$, $\Delta \mathbf{u}_1 = \mathbf{u}_1 - \mathbf{u}_0$, $\Delta \mathbf{u}_2 = \mathbf{u}_2 - \mathbf{u}_0$. We want $\mathbf{p}_u, \mathbf{p}_v$ solving the $2 \times 2$ linear system coming from

$$\begin{align*} \Delta \mathbf{p}_1 &= \mathbf{p}_u \Delta \mathbf{u}_{1_u} + \mathbf{p}_v \Delta \mathbf{u}_{1_v} \\ \Delta \mathbf{p}_2 &= \mathbf{p}_u \Delta \mathbf{u}_{2_u} + \mathbf{p}_v \Delta \mathbf{u}_{2_v} \end{align*}$$

Equivalently,

$$\begin{bmatrix}\Delta \mathbf{p}_1 & \Delta \mathbf{p}_2\end{bmatrix} = \begin{bmatrix}\mathbf{p}_u & \mathbf{p}_v\end{bmatrix}\begin{bmatrix}\Delta \mathbf{u}_{1_u} & \Delta \mathbf{u}_{2_u} \\ \Delta \mathbf{u}_{1_v} & \Delta \mathbf{u}_{2_v}\end{bmatrix}$$

Solving for $\begin{bmatrix}\mathbf{p}_u & \mathbf{p}_v\end{bmatrix}$ by inverting the $2 \times 2$ UV matrix, the determinant is

$$\Delta = \Delta \mathbf{u}_{1_u} \Delta \mathbf{u}_{2_v} - \Delta \mathbf{u}_{1_v} \Delta \mathbf{u}_{2_u}$$

If $\Delta \neq 0$, the inverse yields the closed form

$$\mathbf{p}_u = \frac{\mathbf{u}_{2_v}\mathbf{p}_1 - \mathbf{u}_{1_v}\mathbf{p}_2}{\Delta}, \quad \mathbf{p}_v = \frac{-\mathbf{u}_{2_u}\mathbf{p}_1 - \mathbf{u}_{1_u}\mathbf{p}_2}{\Delta}$$

If $|\Delta| \leq \varepsilon$, the UV mapping is degenerate (no well-defined mapping), so we fall back. Next we Gram-Schmidt the tangent. We set Vector3f n = its.shFrame.n;, then to compute $t$, we project $\mathbf{p}_u$ on to normal $\mathbf{n}$ by $\text{proj}_\mathbf{n}(\mathbf{p}_u) = (\mathbf{n} \cdot \mathbf{p}_u)\mathbf{n}$. Then remove the normal component

$$\tilde{\mathbf{t}} = \mathbf{p}_u - \text{proj}_\mathbf{n}(\mathbf{p}_u) = \mathbf{p}_u - (\mathbf{n} \cdot \mathbf{p}_u)\mathbf{n}$$

and normalize it. Then compute bitangent as cross product: $\mathbf{b} = \mathbf{n} \times \mathbf{t}$. Then we do a linear change-of-basis from tangent-space coordinates $\mathbf{n}_s$ to world coordinates using the TBN matrix $\begin{bmatrix}\mathbf{t} & \mathbf{b} & \mathbf{n}\end{bmatrix}$:

$$\mathbf{n}_\text{world} = n_x \mathbf{t} + n_y \mathbf{b} + n_z \mathbf{n}$$

and normalize it. Then, we replace the shading normal used for shading with the perturbed normal $\mathbf{n}_\text{world}$. The new shading frame is constructed from $\mathbf{n}_\text{world}$.

For the fallback branch (degenerate UV) or if there is no UV, we cannot reliably compute $\mathbf{p}_u, \mathbf{p}_v$ from UV, so we use the preexisting tangent frame.

Validation ↑

The materials used for validation are downloaded from ambientCG:

Comparison: Bark and Ground Material With and Without Normal Maps

Light Transport Algorithms ↑

Final Gather for Photon Mapping ↑

Final Gather for Photon Mapping refines global illumination by gathering photon data to compute indirect light.

In standard photon mapping, we estimate irradiance by averaging photons within a radius $r$. This is inherently biased because we are averaging over an area, not evaluating at a point. Sharp features get smeared. The larger the radius, the more blur, and the smaller the radius, the more noise. Final gather hides the bias behind an integration.

Description ↑

Files created or modified:

• src/photonmapper.cpp

Validation scenes:

• scenes/nori-beyond/final_gather/cbox_0.xml
• scenes/nori-beyond/final_gather/cbox_64_10M.xml

The basic photon mapping logic has two phases: preprocessing and rendering. The preprocessing phase follows the steps:

1. Emit photons from light sources until the desired number of photons stored in map num_photons reaches the target number of photons m_photonCount. Randomly select an emitter, sample a photon ray and power using samplePhoton().
2. Trace photons through the scene :
   • When a photon hits a diffuse surface, store it in the photon map (position, incident direction, power)
   • Use Russian roulette for path termination
   • Sample BSDF to bounce the photon, attenuating its power by $f_\text{r}(\mathbf{x}, \boldsymbol{\omega}_\text{i}, \boldsymbol{\omega}_\text{r}) \cos \theta_\text{i} / p(\boldsymbol{\omega}_i)$
   • Account for radiance scaling when passing through dielectric surfaces with $\eta^2$ scaling
3. Build KD-tree for efficient photon lookups

The rendering phase follows the steps:

1. Trace camera rays through the scene
2. If the ray hits an emitter, add emitted radiance $L_\text{e}$.
3. At the first diffuse surface hit,
   • Search photon map within m_photonRadius
   • For each found photon $i$, evaluate BSDF (with the current ray direction $\boldsymbol{\omega}_\text{o}$ and the photon direction $\boldsymbol{\omega}_\text{i}$) and accumulate $L_\text{r}$ with $f_\text{r}(\mathbf{x}, \boldsymbol{\omega}_\text{i}, \boldsymbol{\omega}_\text{o}) \Phi_i$.
   • Get density estimation $L_\text{r} \leftarrow L_\text{r} / (\pi r^2 M)$ where $M$ is the number of emitted photons.
   • Aggregate $L_\text{i}$ by $L_\text{r}$ multiplied by the current throughput
4. Use Russian roulette for path termination
5. At specular surfaces, continue tracing via BSDF sampling.

PhotonMapper::Li()

Instead of looking up photons at the point of query $\mathbf{x}$, we do one extra bounce. We sample m_numGatherRays gather rays from $\mathbf{x}$ with its BSDF. We trace each gather ray through the scene with max maxSpecularBounces specular bounces. At the first diffuse hit, we perform the photon lookup and density estimation. The radiance $\mathbf{x}$, $L_\text{r}(\mathbf{x}, \boldsymbol{\omega}_\text{o})$ from final gather is then given by

$$L_\text{r}(\mathbf{x}, \boldsymbol{\omega}_\text{o}) = \frac{1}{\texttt{m\_numGatherRays}}\sum_{k = 1}^\texttt{m\_numGatherRays} f_\text{r}(\mathbf{x}, \boldsymbol{\omega}_{\text{i}, k}, \boldsymbol{\omega}_\text{o})L_\text{i}(r(\mathbf{x}, \boldsymbol{\omega}_{\text{i}, k}), -\boldsymbol{\omega}_{\text{i}, k})$$

A caveat about final gather is that we need to explicitly compute direct lighting with a shadow ray test. Consider direct photon density estimation: Without final gather, photons from the light hits surfaces and get stored. A sphere blocks photons from reaching the floor behind it. Fewer photons are then stored in that shadowed region, causing it to be darker and creating shadows. With final gather, gather rays are shot from the shading point in all directions, where we look up photons. Even if the shading point is in shadow, the gather ray can reach lit surfaces, causing the shadow to disappear. Hence the radiance at $\mathbf{x}$, $L_\text{i}(\mathbf{x}, \boldsymbol{\omega}_\text{o})$ is given by $L_\text{i}(\mathbf{x}, \boldsymbol{\omega}_\text{o}) = L_\text{e} + qL_\text{r}(\mathbf{x}, \boldsymbol{\omega}_\text{o})$ where $q$ is the throughput.

The current code now takes in an argument numGatherRays as follows:

<integrator type="photonmapper">
		<float name="photonRadius" value="0.05"/>
		<integer name="photonCount" value="10000000"/>
        <float name="numGatherRays" value="64"/>
</integrator>

If numGatherRays is unspecified, it is automatically set to 0, and the behaviour is the same as original version of photon mapping.

Validation ↑

It can be observed that the caustic effect gets diluted or lost. This is because caustic photons (light paths like $LS+D$) are stored in the global photon map. Without final gather, we query photons directly at the shading point, and caustics are visible. With final gather, the gather rays sample the BSDF and do photon lookups at the endpoints of those rays. However, caustics are highly localized and the probability that a random gather ray lands exactly in a caustic region is very low, and thus most gather rays missthe caustic entirely. So the caustic effect gets diluted.

Dedicated Caustics Photon Map ↑

A dedicated caustics map separates photons into two categories: Global photon map and caustics photon map. The caustics photon map stores photons that passed through at least one specular surface before hitting a diffuse surface.

Description ↑

Files created or modified:

• src/photonmapper.cpp

Validation scenes:

• scenes/nori-beyond/caustics_map/cbox_64_10M.xml

New member variables m_causticPhotonCount, m_causticRadius, m_emittedCausticPhotons have been added. In general, typical photon radius in the caustics map is smaller than that in the global map, and fewer photons are needed.

PhotonMapper::preprocess()

The preprocessing phase now needs to build two photon maps. photonCount photons are to be stored in the global photon map and stored at diffuse surfaces. If caustics photon map is enabled, the photons are stored only if they did not just hit a specular surface one bounce ago. Otherwise if caustics photon map is disabled, they are stored at every diffuse hit. causticPhotonCount photons are to be stored in the caustics photon map. We track a flag hasHitSpecular, initialized as false, as the photon bounces in the scene. Once the photon has hit a specular surface, the flag is set to true. In the next diffuse hit, the photon is stored in the caustics map and we stop tracing this photon.

PhotonMapper::Li()

At a diffuse surface, we estimate the photon density from the caustics map without final gather and obtain $L_\text{caustic}$. We also estimate the photon density with the global photon map with or without final gather (depending on user specification) and obtain $L_\text{global}$. We then aggregate the two radiances.

The current code now takes in an argument causticPhotonCount and causticRadius as follows:

<integrator type="photonmapper">
    <integer name="photonCount" value="10000000"/>
    <float name="photonRadius" value="0.05"/>
    <integer name="numGatherRays" value="64"/>

    <integer name="causticPhotonCount" value="2000000"/>
    <float name="causticRadius" value="0.02"/>
</integrator>

If causticPhotonCount or causticRadius are unspecified, they are automatically set to 0, and the behaviour is the same as original version of photon mapping.

Validation ↑

Resampled Importance Sampling ↑

Resampled Importance Sampling (RIS) focuses computation on areas of higher light contribution.

In the current path tracing with multiple importance sampling (MIS), we sample one light uniformly from the light set, and one BSDF direction, and combine the two with MIS weights. This works well, but uniform light selection is bad when many lights exist, and one light sample per bounce is often noisy. In RIS, we generate $M$ candidate lights, assign each a weight, and resample one candidate. When $M = 1$, RIS reduces to standard emitter sampling.

Description ↑

Files created or modified:

• src/path_ris.cpp

Validation scenes:

• scenes/nori-beyond/path_ris/sponza.xml

The general RIS framework from Talbot et al. estimates an integral $I = \int_\Omega f(x) \, d\mu(x)$ using two PDFs:

• Source PDF $p(x)$: easy to sample from (the proposal distribution)
• Target PDF $g(x)$: closer to $f(x)$, possibly unnormalized or difficult to sample directly

The RIS procedure is:

1. Draw $M$ candidate samples $X_1, \ldots, X_M \sim p(x)$
2. Compute resampling weights $w_j = \frac{g(X_j)}{p(X_j)}$
3. Resample one candidate $Y$ from $\{X_j\}$ with probability $P(Y = X_j) = \frac{w_j}{\sum_k w_k}$
4. The unbiased RIS estimator is:

$$\hat{I}_\text{RIS} = \underbrace{\frac{f(Y)}{g(Y)}}_\text{evaluated at resampled $Y$} \cdot \frac{1}{M} \sum_{j=1}^M w_j$$

When $M = 1$, RIS reduces to standard importance sampling with PDF $p$.

Our implementation uses a simplified parameterization where we set $g = f$ (the target equals the integrand). This is a valid special case of RIS.

For direct lighting at shading point $\mathbf{x}$, the integrand is:

$$f(\mathbf{y}) = f_\text{r}(\mathbf{x}, \boldsymbol{\omega}_\text{i}, \boldsymbol{\omega}_\text{o}) \, L_\text{e}(\mathbf{y}, -\boldsymbol{\omega}_\text{i}) \, \cos\theta_\text{i} \, G(\mathbf{x}, \mathbf{y}) \, V(\mathbf{x}, \mathbf{y})$$

where $\mathbf{y}$ is a point on a light source, $\boldsymbol{\omega}_\text{i} = \frac{\mathbf{y} - \mathbf{x}}{\|\mathbf{y} - \mathbf{x}\|}$, $G$ is the geometry term, and $V$ is visibility.

Our source (proposal) PDF is uniform light selection combined with area sampling on each light:

$$p(\mathbf{y}) = \frac{1}{N_\text{lights}} \cdot p_A(\mathbf{y} \mid l)$$

with $p_A(\mathbf{y} \mid l) = \frac{\cos\theta_\text{o}}{\|\mathbf{x} - \mathbf{y}\|^2} p_\Omega(\boldsymbol{\omega}_\text{i})$.

PathRISIntegrator::Li()

At shading point $\mathbf{x}$:

1. Sample $M$ candidates. For $i = 1, \ldots, M$, sample light $l_i$ uniformly such that $p(l_i) = \frac{1}{N_\text{lights}}$, and sample point $\mathbf{y}^i$ on light $l_i$: $\mathbf{y}^i \sim p_A(\mathbf{y} \mid l_i)$

2. For each candidate, compute the contribution (ignoring visibility for now):

$$c_i = \frac{f_\text{r} \, L_\text{e} \, \cos\theta_\text{i}}{p_\Omega(\boldsymbol{\omega}_\text{i})}$$

Note that $c_i$ is the integrand $f$ divided by the solid angle PDF only (not the full proposal PDF). The $\frac{1}{N_\text{lights}}$ factor is handled separately. Also compute the scalar resampling weight $w_i = \text{luminance}(c_i)$.

3. Select candidate $k$ with probability $P(k) = \frac{w_k}{\sum_{j=1}^M w_j}$

4. If the selected candidate passes the visibility test, compute the final estimate. When $g = f$ (our case), $\frac{f(Y)}{g(Y)} = 1$, so the estimator simplifies to $\frac{1}{M}\sum_j w_j$. However, computing the full sum requires evaluating all $M$ candidates. In practice, we select candidate $k$ with probability $P(k) = \frac{w_k}{\sum w_j}$, then output $\frac{1}{P(k)} \cdot \frac{1}{M} \cdot \frac{f(X_k)}{g(X_k)}$. By importance sampling, this has the same expectation as the full sum but only requires evaluating one candidate (e.g., one visibility test instead of $M$). Using the this single-sample RIS variant with $g = f$, the estimator is:

$$\frac{1}{P(k)} \cdot \frac{1}{M} \cdot \underbrace{\frac{f(X_k)}{g(X_k)}}_{= 1 \text{ since } g = f} = \frac{\sum w_j}{w_k} \cdot \frac{1}{M}$$

But this assumes $w_j = \frac{f(X_j)}{p(X_j)}$ with the full proposal PDF $p = \frac{p_\Omega}{N_\text{lights}}$. In our implementation, we instead use $w_j = \text{luminance}(c_j)$ where $c_j = \frac{f_j}{p_\Omega}$ (missing the $N_\text{lights}$ factor). To correct for this, we multiply by $c_k \cdot N_\text{lights}$:

$$\hat{L}_\text{RIS} = \frac{1}{M} \cdot \frac{\sum_{j=1}^M w_j}{w_k} \cdot \underbrace{c_k \cdot N_\text{lights}}_{= f_k / p(X_k)}$$

The term $c_k \cdot N_\text{lights} = \frac{f_k}{p_\Omega} \cdot N_\text{lights} = \frac{f_k}{p_\Omega / N_\text{lights}} = \frac{f_k}{p(X_k)}$ gives us the correct importance sampling ratio for the full proposal PDF.

When $M = 1$, $\frac{\sum w_j}{w_k} = 1$, recovering standard emitter sampling: $\hat{L} = c_1 \cdot N_\text{lights} = \frac{f}{p}$.

For the emitter sampling branch (RIS), we apply a MIS weight:

$$w_\text{em} = \frac{p_\text{em}}{p_\text{em} + p_\text{mat}}$$

where $p_\text{em} = \frac{p_\Omega(\boldsymbol{\omega}_\text{i})}{N_\text{lights}}$ is the full emitter sampling PDF (solid angle PDF divided by number of lights) and $p_\text{mat}$ is the BSDF sampling PDF. The final RIS contribution becomes:

$$\hat{L}_\text{RIS} = w_\text{em} \cdot \frac{1}{M} \cdot \frac{\sum w_j}{w_k} \cdot c_k \cdot N_\text{lights}$$

For the BSDF sampling branch, when a BSDF-sampled direction hits an emitter, we apply the complementary MIS weight:

$$w_\text{mat} = \frac{p_\text{mat}}{p_\text{em} + p_\text{mat}}$$

and add:

$$\hat{L}_\text{BSDF} = w_\text{mat} \cdot \frac{f_\text{r} \cos\theta}{p_\text{mat}} \cdot L_\text{e}$$

For delta BSDFs (mirrors, glass), $p_\text{mat}$ is a Dirac delta and emitter sampling cannot find specular paths, so we use $w_\text{mat} = 1$ for BSDF-sampled paths through specular surfaces.

Validation ↑

This Sponza scene features multiple area lights distributed throughout the space, creating complex lighting conditions where uniform light selection becomes inefficient. Many shading points receive significant contribution from only a subset of the lights (e.g., nearby lights or those with unobstructed paths), while distant or occluded lights contribute little.

RIS excels in this scenario by resampling candidates based on their actual contribution (BSDF $\times$ emission $\times$ geometry), effectively concentrating samples on the lights that matter most for each shading point. This leads to lower variance compared to MIS with uniform light selection at the same sample count.

Participating Media ↑

Phase Functions ↑

The phase function $f_\text{p}(\mathbf{x}, \boldsymbol{\omega}_\text{i}, \boldsymbol{\omega})$ describes the angular distribution of light as it is scattered at a point within a medium. It is the volumetric analog of the BRDF for surfaces. Unlike a BRDF, which is normalized by a cosine factor and can integrate to values less than 1 (due to absorption), a phase function represents a pure probability density function (PDF) over the sphere of directions. Therefore, it must integrate to unity:

$$\int_{S^2} f_\text{p}(\mathbf{x}, \boldsymbol{\omega}_\text{i}, \boldsymbol{\omega}) \mathrm{d} \boldsymbol{\omega}_\text{i} = 1$$

Description ↑

The phase functions implemented are the isotropic phase function and the henyey-greenstein phase function.

In isotropic scattering, light is scattered equally in all directions. This is the volumetric equivalent of a perfectly diffuse Lambertian surface. Since there are $4\pi$ steradians in a unit sphere,

$$f_\text{p}(\boldsymbol{\omega}_\text{i}, \boldsymbol{\omega}_\text{o}) = \frac{1}{4\pi}$$

In many real-world media (like clouds or smoke), scattering is highly anisotropic, meaning it favors certain directions. Henyey-Greenstein is a popular empirical model used to simulate this behavior:

$$f_{\text{p}, \text{HG}}(\theta) = \frac{1}{4\pi} \frac{1 - g^2}{(1 + g^2 - 2g \cos \theta)^{3/2}}$$

where $\theta$ is the angle between the incoming and outgoing directions. The parameter $g$ (the asymmetry parameter or mean cosine) defines the "shape" of the scattering:

• $g > 0$: Forward scattering (light continues mostly in the same direction).
• $g < 0$: Backscattering (light is reflected back toward the source).
• $g = 0$: Isotropic scattering (mathematically equivalent to the isotropic model).

The value of $g$ is defined as the average cosine of the scattering angle:

$$g = \int_{S^2} f_\text{p}(\mathbf{x}, \boldsymbol{\omega}_\text{i}, \boldsymbol{\omega}) \cos \theta \, \mathrm{d} \boldsymbol{\omega}_\text{i}$$

where $\cos \theta = -\boldsymbol{\omega} \cdot \boldsymbol{\omega}_\text{i}$. For example, if $g=0.5$, the "average" photon is deflected only $60^\circ$ ($\cos 60^\circ = 0.5$) from its original path.

Files created or modified:

• include/nori/phasefunction.h
• src/isotropic.cpp
• src/henyey_greenstein.cpp

Validation scenes:

include/nori/phasefunction.h

This file mirrors the BSDF interface for consistency. It includes a PhaseFunctionQueryRecord, similar to the BSDFQueryRecord, and supports two constructors of a PhaseFunctionQueryRecord -- one that takes in the incident direction $\boldsymbol{\omega}_\text{i}$ for sampling the outgoing direction $\boldsymbol{\omega}_\text{o}$, and one that takes in both $\boldsymbol{\omega}_\text{i}$ and $\boldsymbol{\omega}_\text{o}$ for evaluation.

The PhaseFunction object is the superclass of all phase functions and has the methods sample and eval. Unlike BSDF which also has the pdf method, phase functions are PDFs themselves, hence eval will return the phase function value (i.e. the PDF). The sample method in PhaseFunction also returns the phase function value for convenience (instead of 1.0, which would be case in the BSDF implementation that returns eval() / pdf()).

Isotropic::sample(...)

Given a sample $\boldsymbol{\xi} \in [0, 1]^2$, we sample a direction uniformly on a sphere with Warp::squareToUniformSphere and assign this direction to pRec.wo. The method returns the phase function value is $1 / (4\pi)$.

Isotropic::eval(...)

The phase function value is simply $1 / (4\pi)$ for any given pairs of $\boldsymbol{\omega}_\text{i}$ and $\boldsymbol{\omega}_\text{o}$.

HenyeyGreensteinPhaseFunction::sample(...)

The azimuthal angle is simply $\phi = 2\pi \xi_2$ since Henyey-Greenstein does not depend on it. For the polar angle, we need to find $\cos\theta$ by inverting the CDF. The marginal PDF for $\cos\theta$ is obtained by integrating over $\phi$:

$$p(\cos\theta) = \frac{1 - g^2}{2(1 + g^2 - 2g\cos\theta)^{3/2}}.$$

The CDF is obtained by integrating from $-1$ to $\cos\theta$. Let $u = 1 + g^2 - 2g\mu$, so $\mathrm{d}u = -2g \, \mathrm{d}\mu$:

$$\begin{align*} P(\cos\theta) &= \int_{-1}^{\cos\theta} \frac{1 - g^2}{2(1 + g^2 - 2g\mu)^{3/2}} \, \mathrm{d}\mu \\ &= \frac{1 - g^2}{2} \int_{(1+g)^2}^{1 + g^2 - 2g\cos\theta} u^{-3/2} \cdot \frac{\mathrm{d}u}{-2g} \\ &= \frac{1 - g^2}{-4g} \left[ -2u^{-1/2} \right]_{(1+g)^2}^{1 + g^2 - 2g\cos\theta} \\ &= \frac{1 - g^2}{2g} \left( \frac{1}{\sqrt{1 + g^2 - 2g\cos\theta}} - \frac{1}{1+g} \right). \end{align*}$$

Setting $P(\cos\theta) = \xi_1$ and solving for $\cos\theta$:

$$\begin{align*} \xi_1 &= \frac{1 - g^2}{2g} \left( \frac{1}{\sqrt{1 + g^2 - 2g\cos\theta}} - \frac{1}{1+g} \right) \\ \frac{2g\xi_1}{1 - g^2} + \frac{1}{1+g} &= \frac{1}{\sqrt{1 + g^2 - 2g\cos\theta}} \\ \frac{2g\xi_1 + (1-g)}{1 - g^2} &= \frac{1}{\sqrt{1 + g^2 - 2g\cos\theta}} \\ \frac{1 - g + 2g\xi_1}{(1-g)(1+g)} &= \frac{1}{\sqrt{1 + g^2 - 2g\cos\theta}} \\ \sqrt{1 + g^2 - 2g\cos\theta} &= \frac{(1-g)(1+g)}{1 - g + 2g\xi_1} = \frac{1 - g^2}{1 + g - 2g(1 - \xi_1)}. \end{align*}$$

Squaring both sides and solving for $\cos\theta$:

$$\cos\theta = \frac{1}{2g}\left(1 + g^2 - \left(\frac{1 - g^2}{1 + g - 2g\xi_1}\right)^2\right).$$

For the special case $g \approx 0$ (isotropic), we use uniform sphere sampling: $\cos\theta = 1 - 2\xi_1$.

Once $\cos\theta$ and $\phi$ are computed, we build a coordinate frame around the forward scattering direction $-\boldsymbol{\omega}_\text{i}$ and convert from spherical to Cartesian coordinates to obtain $\boldsymbol{\omega}_\text{o}$.

HenyeyGreensteinPhaseFunction::eval(...)

The method simply computes $\cos \theta$ for the given pair of $\boldsymbol{\omega}_\text{i}$ and $\boldsymbol{\omega}_\text{o}$, and returns

$$f_{\text{p}, \text{HG}}(\theta) = \frac{1}{4\pi} \frac{1 - g^2}{(1 + g^2 - 2g \cos \theta)^{3/2}}.$$

Validation ↑

Transmittance in Homogeneous Media ↑

This feature enables coloured dielectric materials and fog volumes using Beer's law transmittance for homogeneous participating media.

Description ↑

By Beer's law, in a homogeneous participating medium (where the absorption coefficient $\boldsymbol{\sigma}_\text{a}$ is constant), the transmittance of light is given by

$$T = \exp\{-\boldsymbol{\sigma}_\text{a} t\}$$

where $t$ is the distance travelled.

Files created or modified:

• include/nori/bsdf.h
• src/null.cpp
• src/dielectric.cpp
• src/path_mis.cpp

Validation scenes:

• scenes/nori-beyond/homogeneous_media/cbox_colored_glass.xml
• scenes/nori-beyond/homogeneous_media/cbox_colored_glass_mitsuba.xml

include/nori/bsdf.h

The BSDF interface is extended to include absorption with default no-absorption behaviour. New virtual methods getAbsorptionCoefficient() (returns Color3f(0.0f) by default) and hasAbsorption() (returns false by default) are added. A method isNull() (returns false by default) is also added to distinguish null/transparent surfaces from real surfaces during ray marching.

src/null.cpp

This is a null/transparent BSDF that lets rays pass through without any scattering or reflection. It is used to bound fog volumes without any refractive distortion. It overrides hasAbsorption and getAbsorptionCoefficient to expose the volume's absorption coefficient, and isNull to return true. The null BSDF has discrete measure and returns zero for eval and pdf. In the sample method, it sets $\boldsymbol{\omega}_\text{o}$ to $-\boldsymbol{\omega}_\text{i}$ (ray passes straight through) and $\eta$ to 1 (no index of refraction change). It can be used in the scene like this:

<mesh type="obj">
    <string name="filename" value="fog_box.obj"/>
    <bsdf type="null">
        <color name="sigmaA" value="0.5, 0.5, 0.5"/>
    </bsdf>
</mesh>

src/dielectric.cpp

The dielectric BSDF is extended with a new Color3f member variable m_sigmaA, which stores the absorption coefficient $\boldsymbol{\sigma}_\text{a}$ per RGB channel. The methods getAbsorptionCoefficient and hasAbsorption are overridden. The absorption itself is not applied by the BSDF; it is handled by the integrator via Beer's law on each ray segment inside the medium.

PathMISIntegrator::Li(...)

The integrator uses a medium stack (std::vector``<``const BSDF *>) to track which absorbing volumes the ray is currently inside. This stack-based approach correctly handles arbitrarily nested media (e.g. a coloured glass sphere inside a fog cube). The innermost medium on the stack determines the active absorption coefficient.

The core procedure is null-surface marching: given a ray, we step through all null (transparent) boundaries until we reach a real surface, applying Beer's law and updating the medium stack at each step. This procedure is used in three contexts:

1. Initial ray marching (camera ray to first real surface)
2. Shadow rays (emitter sampling, using a separate copy of the medium stack)
3. BSDF-sampled rays (path continuation after bouncing off a real surface)

The null-surface marching procedure works as follows:

1. Trace the ray to the nearest surface intersection.
2. Look up the current absorption coefficient $\boldsymbol{\sigma}_\text{a}$ from the top of the medium stack (zero if the stack is empty).
3. Apply Beer's law: multiply the throughput by $T = \exp\{-\boldsymbol{\sigma}_\text{a} \cdot t\}$, where $t$ is the distance to the intersection.
4. If the intersected surface is not null (i.e. it is a real surface like a diffuse wall or dielectric), stop marching and return this intersection.
5. Otherwise, the surface is null (a medium boundary). Update the medium stack:
   • If the ray is entering the volume (ray direction dot geometric normal $< 0$), push the null surface's BSDF onto the stack.
   • If the ray is exiting the volume (ray direction dot geometric normal $> 0$), find and remove this BSDF from the stack.
6. Advance the ray origin to the intersection point and go back to step 1.

For shadow rays, the procedure is identical except that it uses a separate copy of the medium stack (so it does not disturb the main path state), accumulates transmittance into a separate variable, and terminates early if a non-null surface is hit (the light is occluded).

For dielectric surfaces, the medium stack is updated when the BSDF samples a transmission (refraction) event. This is detected by checking whether $\boldsymbol{\omega}_\text{i}$ and $\boldsymbol{\omega}_\text{o}$ are on opposite sides of the geometric normal. If so, the ray has crossed the interface and the dielectric's absorption medium is pushed onto or removed from the stack.

The medium stack uses id-based removal: when exiting a volume, the stack is searched from the back (innermost) to find and remove the matching BSDF pointer. This is necessary because overlapping volumes do not always nest cleanly. In the validation scene below, the sphere straddles the fog cube boundary (it pokes above the fog). Consider a ray traveling downward through the top of the sphere into the fog:

1. Ray enters the sphere (above the fog) $\to$ stack: [sphere]
2. Ray enters the fog cube (partway through the sphere) $\to$ stack: [sphere, fog]
3. Ray exits the sphere (still inside the fog) $\to$ need to remove sphere, but fog is on top

At step 3, a naive stack pop would remove fog instead of sphere, which is wrong. Searching from the back by BSDF identity correctly finds and removes sphere from below fog, leaving the stack as [fog].

Validation ↑

The validation scene places a blue glass sphere (null BSDF, $\boldsymbol{\sigma}_\text{a} = (3.0, 3.0, 0.1)$) and an amber glass sphere (dielectric BSDF, $\boldsymbol{\sigma}_\text{a} = (0.5, 1.5, 4.0)$) inside a fog cube (null BSDF, $\boldsymbol{\sigma}_\text{a} = (0.5, 0.5, 0.5)$) in a Cornell box. The fog cube extends through the bottom part of the room, while the spheres poke above the fog, so the fog boundary crosses each sphere near its equator.

The Nori and Mitsuba renders differ due to a fundamental difference in how the two renderers track media. Mitsuba 3 uses a single medium pointer per ray: at each surface crossing, the active medium is replaced with the surface's interior or exterior medium. This produces two visible artifacts in the Mitsuba renders:

1. A visible horizontal plane across the spheres at the fog cube boundary. This appears as a brightness discontinuity where the fog cube's top face intersects the spheres. Because Mitsuba's medium pointer switches discretely at the null boundary, there is an abrupt transition between "fog active" (below the plane) and "no fog" (above the plane). In Nori, the stack-based approach produces a smooth transition because entering/exiting the fog and sphere are tracked independently.

2. A bright ring at the sphere-fog boundary in the original Mitsuba render (center image). When a ray exits the sphere (still geometrically inside the fog cube), Mitsuba sets the medium to the sphere's exterior medium. If no exterior medium is specified on the sphere, Mitsuba treats the region as vacuum, losing the fog absorption for the remainder of the segment. Rays grazing the sphere edge barely enter the sphere (little sphere absorption) and also miss the fog behind it (medium reset to vacuum on exit), producing the bright ring.

Nori uses a stack-based approach: entering a volume pushes its BSDF onto the stack, and exiting pops it. When a ray exits the sphere, the fog is still on the stack, so fog absorption continues to be applied correctly. This produces a smoother, physically correct result without the horizontal plane artifact or the bright ring (left image).

To verify this, I re-rendered the Mitsuba scene with <ref id="fog_medium" name="exterior"/> added to both spheres, which tells Mitsuba to restore the fog medium on sphere exit (right image). This largely eliminates the bright ring artifact and brings Mitsuba closer to Nori's output. However, Mitsuba's single-pointer approach still cannot fully represent the overlapping geometry: the sphere's exterior medium is set to fog everywhere, including the upper hemisphere that is above the fog cube. This means Mitsuba now over-attenuates above the fog boundary while correctly attenuating below it. Nori's stack naturally handles this because it tracks all volume boundaries the ray has crossed, regardless of the order.