CS295 Realistic Image Synthesis
Programming Assignment 1: Monte Carlo Sampling and Direct Illumination

In this exercise you will generate sample points on various domains: the plane, disks, spheres, and hemispheres. The base code has been extended with an interactive visualization and testing tool to make working with point sets as intuitive as possible.

After compiling the provided base code, you should see an executable named warptest. Run this executable to launch the interactive warping tool, which allows you to visualize the behavior of different warping functions given a range of input point sets (independent, grid, and stratified). Up to now, we only discussed uniform random variables which correspond to the "independent" type, and you need not concern yourself with the others for now.

Part 1 is split into several subsections; in each case, you are asked to implement a distribution function and a matching sample warping scheme It is crucial that both are consistent with respect to each other (i.e. that warped samples have exactly the distribution described by the density function). Significant errors can arise if inconsistent warpings are used for Monte Carlo integration. The warptest tool provided by us implements a \(\chi^2\) test to ensure that this consistency requirement is indeed satisfied.

Part 1.1. Sample Warping

Implement the missing functions in class Warp found in src/warp.cpp. This class consists of various warp methods that take as input a 2D point \((s, t) \in [0, 1) \times [0, 1) \) and return the warped 2D (or 3D) point in the new domain. Each method is accompanied by another method that returns the probability density with which a sample was picked. Our default implementations all throw an exception, which produces an error message in the graphical user interface of warptest. The slides on the course website provide a number of useful recipes for warping samples and computing the densities, and the PBRT textbook also contains considerable information on this topic that you should feel free to use.

• Warp::squareToTent and Warp::squareToTentPdf

  Implement a method that transforms uniformly distributed 2D points on the unit square into the 2D "tent" distribution, which has the following form: \[ p(x, y)=p_1(x)\,p_1(y)\quad\text{and}\quad p_1(t) = \begin{cases} 1-|x|, & -1\le x\le 1\\ 0,&\text{otherwise}\\ \end{cases} \]

  Note that this distribution is composed of two independent 1D distributions, which makes this task considerably easier. Follow the "recipe" discussed in class:
  1. Compute the CDF \(P_1(t)\) of the 1D distribution \(p_1(t)\)
  2. Derive the inverse \(P_1^{-1}(t)\)
  3. Map a random variable \(\xi\) through the inverse \(P_1^{-1}(t)\) from the previous step
  Show the details of these steps in your report (either using TeX, or by taking a photograph of the derivation and embedding the image)

• Warp::squareToUniformDisk and Warp::squareToUniformDiskPdf

  Implement a method that transforms uniformly distributed 2D points on the unit square into uniformly distributed points on a planar disk with radius 1 centered at the origin. Next, implement a probability density function that matches your warping scheme.

• Warp::squareToUniformSphere and Warp::squareToUniformSpherePdf

  Implement a method that transforms uniformly distributed 2D points on the unit square into uniformly distributed points on the unit sphere centered at the origin. Implement a matching probability density function.

• Warp::squareToUniformHemisphere and Warp::squareToUniformHemispherePdf

  Implement a method that transforms uniformly distributed 2D points on the unit square into uniformly distributed points on the unit hemisphere centered at the origin and oriented in direction \((0, 0, 1)\). Add a matching probability density function.

• Warp::squareToCosineHemisphere and Warp::squareToCosineHemispherePdf

  Transform your 2D point to a point distributed on the unit hemisphere with a cosine density function \[ p(\theta)=\frac{\cos\theta}{\pi}, \] where \(\theta\) is the angle between a point on the hemisphere and the north pole.

Part 1.2. Validation

Pass the \(\chi^2\) test for each one of the above sampling techniques and include screen shots in your report.

In this part of the homework, you'll implement two basic rendering algorithms that set the stage for fancier methods investigated later in the course. For now, both of the methods assume that the object is composed of a simple white diffuse material that reflects light uniformly into all directions.

Part 2.1. Point lights

The provided base code includes a scene scenes/pa1/ajax-simple.xml that instantiates a (currently nonexistent) integrator/rendering algorithm named simple, which simulates a single point light source located at a 3D position position, and which emits an amount of energy given by the parameter energy.

<!-- An excerpt from scenes/pa1/ajax-simple.xml: -->
<integrator type="simple">
    <point name="position" value="-20, 40, 20"/>
    <color name="energy" value="3000, 3000, 3000"/>
</integrator>

Your first task will be to create a new Integrator that accepts these parameters (in a similar way as the dummy normals integrator shown at the beginning of Part 0.4). Take a look at the PropertyList class, which should be used to extract the two parameters.

Let \(\mathbf{p}\) and \(\mathbf{\Phi}\) denote the position and energy of the light source, and suppose that \(\mathbf{x}\) is the point being rendered. Then this integrator should compute the quantity

\[ L(\mathbf{x})=\frac{\Phi}{4\pi^2} \frac{\mathrm{max}(0, \cos\theta)}{\|\mathbf{x}-\mathbf{p}\|^2} V(\mathbf{x}\leftrightarrow\mathbf{p}) \]

where \(\theta\) is the angle between the direction from \(\mathbf{x}\) to \(\mathbf{p}\) and the shading surface normal (available in Intersection::shFrame::n) at \(\mathbf{x}\) and

\[ V(\mathbf{x}\leftrightarrow\mathbf{p}):=\begin{cases} 1,&\text{if $\mathbf{x}$ and $\mathbf{p}$ are mutually visible}\\ 0,&\text{otherwise} \end{cases} \]

is the visibility function, which can be implemented using a shadow ray query. Intersecting a shadow ray against the scene is generally cheaper since it suffices to check whether an intersection exists rather than having to find the closest one.

Implement the simple integrator according to this specification and render the scene scenes/pa1/ajax-simple.xml. Include a comparison against the reference image in your report.

Part 2.2. Ambient occlusion

Ambient occlusion is rendering technique which assumes that a (diffuse) surface receives uniform illumination from all directions (similar to the conditions inside a light box), and that visibility is the only effect that matters. Some surface positions will receive less light than others since they are occluded, hence they will look darker. Formally, the quantity computed by ambient occlusion is defined as

\[ L(\mathbf{x})=\int_{\Omega_\mathbf{x}}V(\mathbf{x}, \mathbf{x}+\alpha\omega)\,\frac{\cos\theta}{\pi}\,\mathrm{d}\omega \]

which is an integral over the upper hemisphere \(\Omega_{\mathbf{x}}\) centered at the point \(\mathbf{x}\). The variable \(\theta\) refers to the angle between the direction \(\omega\) and the shading normal at \(\mathbf{x}\). The ad-hoc variable \(\alpha\) adjusts how far-reaching the effects of occlusion are.

Note that this situation—sampling points on the hemisphere with a cosine weight—exactly corresponds to one of the warping functions you implemented in Part 2.1, specifically squareToCosineHemisphere. Use this function to sample a point on the hemisphere and then check for visibility using a shadow ray query. You can assume that occlusion is a global effect (i.e. \(\alpha=\infty\)).

One potential gotcha is that the samples produced by squareToCosineHemisphere lie in the reference hemisphere and need to be oriented according to the surface at \(\mathbf{x}\). Take a look at the Frame class, which is intended to facilitate this.

Implement the ambient occlusion (ao) integrator and render the scene scenes/pa1/ajax-ao.xml. Include a comparison against the reference image in your report.

Part 3.1. Area lights

Our first goal will be to extend Nori so that any geometric object can be turned into a light source known as an area light.

Each triangle of a mesh that is marked as an area light uniformly emits radiance towards all directions above its surface. In Nori's XML description language, area lights are specified using a nested emitter tag of type area. Here is an example:

<scene>
    <!-- Load a OBJ file named "bunny.obj" -->
    <mesh type="obj">
        <string name="filename" value="bunny.obj"/>

        <!-- Turn the mesh into an area light source -->
        <emitter type="area">
            <!-- Assign a uniform radiance of 1 W/m2sr -->
            <color name="radiance" value="1, 1, 1"/>
        </emitter>
    </mesh>

    <!-- ..... -->
</scene>

Currently, Nori won't be able to understand the above snippet since area lights are not yet implemented. To add area lights to Nori, follow these steps:

1. Create a new class AreaLight in a file named src/area.cpp that derives from the Emitter class. Connect it to the scene parser using the NORI_* macros similar to the Integrators you have previously created. Use the constructor's const PropertyList & argument to extract the radiance parameter in the constructor.
   
   

2. The Monte Carlo rendering technique in Part 3.2 requires the ability to sample points that are uniformly distributed on area lights. Currently, none of this functionality exists.

   Begin by familiarizing yourself with the Mesh class to see how vertices, faces and normals are stored. Next, add a method that uniformly samples positions on the surface associated with a specific Mesh instance. The name and precise interface of this method are completely up to you. However, we suggest that it should take a uniformly 2D sample and return:

   1. The sampled position \(\mathbb{p}\) on the surface of the mesh.
   2. The interpolated surface normal \(\mathbb{n}\) at \(\mathbb{p}\) computed from the per-vertex normals. When the mesh does not provide per-vertex normals, compute and return the face normal instead.
   3. The probability density of the sample. This should be the reciprocal of the surface area of the entire mesh.
      

   You may find the DiscretePDF class (declared in include/nori/dpdf.h) useful to implement the sampling step. We suggest that you use this class to build a discrete probability distribution that will allow you to pick a triangle proportional to its surface area. Once a triangle is chosen, you can (uniformly) sample a barycentric coordinate \((\alpha, \beta, 1-\alpha-\beta)\) using the mapping \[ \begin{pmatrix} \alpha\\ \beta \end{pmatrix} \mapsto \begin{pmatrix} 1 - \sqrt{1 - \xi_1}\\ \xi_2\, \sqrt{1 - \xi_1} \end{pmatrix} \] where \(\xi_1\) and \(\xi_2\) are uniform variates.

   The precomputation to build the discrete probability distribution can be performed in the activate() method of the Mesh class, which is automatically invoked by the XML parser.

Part 3.2. Distribution Ray Tracing

In this part you will implement a new direct illumination integrator, which integrates the incident radiance by sampling points on a set of emitters (a.k.a. light sources). Emitters can be fully, partially or not at all visible from a point in your scene, hence you will need to perform Monte Carlo integration to compute the reflected radiance while accounting for visibility.

Recall the Reflection Equation discussed in class, which expresses the reflected radiance due to incident illumination from all directions as an integral over the unit hemisphere \(\Omega_{\mathbf{x}}\) at \(\mathbf{x}\): \[ \newcommand{\vx}{\mathbf{x}} \newcommand{\vc}{\mathbf{c}} \newcommand{\vy}{\mathbf{y}} \newcommand{\vn}{\mathbf{n}} L_r(\vx,\omega_o) = \int_{\Omega_{\vx}} f_r (\vx,\omega_i \leftrightarrow \omega_o)\,L_i (\vx,\omega_i) \cos\theta_i\, \mathrm{d}\omega_i. \] We'll now put together all of the pieces to approximate this integral using Monte Carlo sampling.

Begin by taking a look at the BSDF class in Nori, which is an abstract interface for materials representing the \(f_r\) term in the above equation. Evaluating \(f_r\) entails a call to the BSDF::eval() function, while sampling and probability evaluation are realized using the BSDF::sample() and BSDF::pdf() methods. All methods take a special BRDFQueryRecord as argument, which stores relevant quantities in a convenient data structure.

In this assignment, we will only consider direct illumination, which means that \(L_i(\vx,\omega_i)\) is zero almost everywhere except for rays that happen to hit an area light source. A correct but naïve way of evaluating this integral would be to uniformly sample a direction on the hemisphere and then check if it leads to an intersection with a light source.

However, doing so would be extremely inefficient: light sources generally only occupy a tiny area on the hemisphere, hence most samples would be wasted, causing the algorithm to produce unusably noisy and unconverged images.

We will thus use a better strategy with a higher chance of success: instead of sampling directions on the hemisphere and checking if they hit a light source, we will directly sample points on the light sources and then check if they are visible as seen from \(\vx\). Conceptually, this means that we will integrate over the light source surfaces \(A_e\) instead of the hemisphere \(\Omega_{\mathbf{x}}\): \[ \newcommand{\vr}{\mathbf{r}} L_r^\mathrm{direct}(\vx,\omega_o) = \int_{A_e} f_r (\vx,(\vx\to\vy) \leftrightarrow \omega_o)\,L_e (\vy,\vy\to\vx) \, \mathrm{d} \vy?? \qquad(\text{warning: this is not (yet) correct}) \]

Here \(\mathbf{x}\to\mathbf{y}\) refers to the normalized direction from \(\mathbf{x}\) to \(\mathbf{y}\) (i.e., \(\mathbf{x}\to\mathbf{y} := (\mathbf{y} - \mathbf{x})/\| \mathbf{y} - \mathbf{x} \|\)) and \(L_e(\vx,\omega)\) is the amount of emitted radiance at position \(\vx\) into direction \(\omega\). The integral above motivates the algorithm, but it is not correct: since we changed the integration variable from the solid angle domain to positions, there should be a matching change of variables factor that accounts for this (this is not unlike switching from polar coordinates to a Cartesian coordinate system). In our case, this change of variable factor is known as the geometric term: \[ G(\vx\leftrightarrow\vy) :=V(\vx, \vy)\frac{ |\vn_\vx \cdot(\vx\to\vy)|\,\cdot\, |\vn_\vy \cdot(\vy\to\vx)|}{\|\vx-\vy\|^2} \]

The first term \(V(\vx, \vy)\) is the visibility function, which is \(1\) or \(0\) if the two points are mutually visible or invisible, respectively. The numerator contains the absolute value of two dot products that correspond to the familiar cosine foreshortening factors at both \(\vx\) and \(\vy\). The denominator is the inverse square falloff that we already observed when rendering with point lights. The "\(\leftrightarrow\)" notation expresses that \(G\) is symmetric with respect to its arguments.

Given the geometric term, we can now write down the final form of the reflection equation defined as an integral over surfaces: \[ \newcommand{\vr}{\mathbf{r}} L_r^\mathrm{direct}(\vx,\omega_o) = \int_{A_e} f_r (\vx,(\vx\to\vy) \leftrightarrow \omega_o)\,G(\vx\leftrightarrow\vy)\,L_e (\vy,\vy\to\vx)\, \mathrm{d} \vy \] Note that the cosine factor in the original integral is absorbed by one of the dot products the geometric term. To implement distribution ray tracing in Nori, follow these steps:

1. Create a new integrator src/whitted.cpp (the name will become clear later)
   
   
2. The integrator should begin by finding the first surface interaction \(\vx\) visible along the ray passed to the Li() function. This part works just like in the simple.cpp integrator. What this step does is to apply the ray tracing operator \(\mathrm{RayTrace}(\vc, \omega_c)\) to convert reflected radiance at surfaces into incident radiance at the camera \(\vc\): \[ L_i(\vc,\omega_c)=L_r(\mathrm{RayTrace}(\vc, \omega_c), -\omega_c). \]

3. Given \(\vx\), the distribution ray tracer should then approximate the above integral by sampling a single position \(\vy\in\mathcal{L}\) and then returning the body of the integral, i.e. \[ f_r (\vx,(\vx\to\vy) \leftrightarrow \omega_o)\,G(\vx\leftrightarrow\vy)\,L_e (\vy,\vy\to\vx) \] divided by the probability of the sample \(\vy\) per unit area. However, this will require a few extra pieces of functionality.

   Take a look at the Emitter interface. It is almost completely empty. Clearly, some mechanism for sample generation, evaluation of probabilities, and for returning the emitted radiance is needed. We don't explicitly specify an API that you should use to implement the sampling and evaluation operations for emitters—finding suitable abstractions is part of the exercise. That said, you can look at the BSDF definitions in include/nori/bsdf.h to get a rough idea as to how one might get started with such an interface.

Discuss the design choices and challenges you faced in your report. Render the scene scenes/pa1/motto-diffuse.xml and don't forget to include a comparison against our reference solution: motto-diffuse.exr.