Assignment 3: Phong and Multiple Importance Sampling

Due November 1, 2017 at 11:59pm on myCourses
Worth 20%

---

Overview

In this assignment, you'll implement BRDF importance sampling for the Phong reflectance model, as well as a direct illumination integrator combining both BRDF and emitter importance sampling using multiple importance sampling (MIS).

---

Task 1: Phong Reflectance Model (50 points)

You'll begin by implementing a new BRDF for the Phong reflectance model. In 1975, Phong introduced a shading model that combines diffuse and glossy reflection. Many more advanced/accurate models exist (e.g., using microfacet theory), however implementing Phong's model in a physically-based renderer remains a nice way to become familiar with what it takes to add different BRDFs to a system.

The Phong BRDF is defined as the sum of a diffuse and a glossy component:

$$\begin{eqnarray} f_r(\mathbf{x}, \omega_i, \omega_r) &=& f_{r,d}(\mathbf{x}, \omega_i, \omega_r) + f_{r,s}(\mathbf{x}, \omega_i, \omega_r) \nonumber \\ &=& \rho_d \frac{1}{\pi} + \rho_s \frac{n+2}{2\pi}\cos^n \alpha, \label{phong} \end{eqnarray}$$

where

• $\rho_d$ is the diffuse reflectivity, the fraction of the incoming energy that is reflected diffusely (clamped to $1$ for energy balance),
• $\rho_s$ is the specular reflectivity, the fraction of the incoming energy that is reflected specularly,
• $n$ is the Phong exponent (higher value give more mirror-like specular reflection),
• $\alpha$ is the specular angle between the perfect specular reflective direction and the lighting direction.

We also require $\rho_d + \rho_s \leq 1$ for energy conservation.

Begin by implementing the Phong BRDF in src/bsdfs/phong.cpp. This class has three properties: the diffuse and specular reflectances $\rho_d$ & $\rho_s$, as well as the Phong exponent $n$ describing the “shininess" of the material. You need to implement the three methods: Phong::eval(), Phong::sample() and Phong::pdf(), detailed below.

BRDF Evaluation (20 points)

Implement Phong::eval() using Equation \eqref{phong}. Be sure to return zero when incident illuminatino is queried from the below the horizong/top-hemisphere about a shading point. Have a look at src/bsdfs/diffuse.cpp to see how this is done.

When complete, render the Phong Stanford bunny scenes/hw3/Tests/bunny-phong-emitter.xml scene. For a change, this result will have color! You should have the following output:

BRDF Importance Sampling (30 points)

In many cases, uniform sampling and emitter sampling will generate high-variance results with a specular/glossy BRDF like Phong. We can help reduce this variance by importance sampling the BRDF. Recall that, when writing Monte Carlo estimates of the reflection equation, we can choose to distribute our samples according to any valid PDF. Choosing a PDF proportional to one of the integrand terms corresponds to importance sampling that term. So far we've only sampled the cosine term or the emitted radiance from light sources; we can also sample proportional to the BRDF such that $p(\omega_i) \propto f_r(\mathbf{x}, \omega_i, \omega_o)$. Let's apply this to the Phong model.

Importance sampling a BRDF first requires that we express the desired distribution in a convenient coordinate system. We can then compute the marginal and conditional 1D PDFs, and sample these PDFs using the inversion method. For Phong, diffuse and specular samples are handled as follows:

Diffuse Hemispherical Sampling

A diffuse sample can be obtained from the cosine-weighted distribution, so simply calling Warp::squareToCosineHemisphere() does the trick, here.

Specular Lobe Sampling

It is impossible to sample the complete specular component of a BRDF directly because of the cosine foreshortening factor; we are, however, able to sample according to $\cos^n\alpha$ factor, which yields the PDF:

\begin{equation} p(\omega_i) = \frac{n+2}{2\pi}\cos^n\alpha. \end{equation}

Following the same marginalization and inversion method, we can arrive at an importance sampled direction within the specular lobe as:

\begin{equation} \omega_i = (\alpha, \phi) = (\arccos\xi_1^{\frac{1}{n+1}}, 2\pi \xi_2), \end{equation}

again, expressed in spherical coordinates. We can convert to Cartesian coordinates for the direction vector as:

$$\begin{eqnarray} (x,y,z) &=& (\sin\alpha \cos\phi, \sin\alpha \sin\phi, \cos\alpha) \\ &=& \left(\sqrt{1 - \xi_1^\frac{2}{n+1}} \cos(2\pi\xi_2), \sqrt{1 - \xi_1^\frac{2}{n+1}} \sin(2\pi\xi_2), \xi_1^\frac{1}{n+1}\right). \end{eqnarray}$$

Sampling Phong

All that remains now is to determine whether we should draw a diffuse or specular sample (note: your MC estimator can be split into two estimators, one for each component, alternatively). To do so, flip a coin $\xi \in [0,1)$ and use the following rule:

• $\xi \in [0, \rho_d)$: take a diffuse sample and compute its contribution,
• $\xi \in [\rho_d, \rho_d + \rho_s)$: take a specular sample and compute its contribution,
• $\xi \in [\rho_d + \rho_s, 1)$: no contribution.

To compute this contribution, you will have to rescale your sample depending on its type. One way to do so is to compute the specular sampling weight $\kappa = \rho_s / (\rho_s + \rho_d)$. If you choose to sample specular, do $\xi \gets \xi/\kappa$. If you sample diffuse otherwise, perform $\xi \gets (\xi-\kappa)/(1-\kappa)$. This will allow you to sample proportionally in your implementation.

Implement the described algorithm in Phong::sample(). Modify the Li() method in your direct integrator src/integrators/direct.cpp to handle the case EMeasure::EBSDF and render scenes/hw3/Tests/bunny-phong-brdf.xml as a test scene.

Warp::EWarpType Warp::getWarpType(const EMeasure measure, const std::string& warpType/* = ""*/) {
    EWarpType result = EWarpType::ENone;

    switch (measure)
    {
    case EMeasure::EUnknownMeasure:
        result = getWarpType(warpType); 
    case EMeasure::ESolidAngle:
        result = EWarpType::EUniformCone;
    case EMeasure::EArea:
        result = EWarpType::EUniformSphere;
    case EMeasure::EHemisphere:
    case EMeasure::EBSDF:
        if (warpType == s_uniformHemisphere)
            result = EWarpType::EUniformHemisphere;
        else
            result = EWarpType::ECosineHemisphere;
    }

    return result; 
}

Below is a reference image to compare against.

---

Task 2: Multiple Importance Sampling (50 points)

In this part of the assignment, you will implement multiple importance sampling (MIS). Specifically, you will implement an integrator that combines BRDF and emitter importance sampling to render an image.

Recall that if two sampling distributions $p_f$ and $p_g$ can be used in an estimator for the integral $\int f(x)g(x)\, \mathrm{d}x$, a single Monte Carlo estimator that uses both distributions simultaneously is given by the MIS formulation as

\begin{equation} \int f(x)g(x)\, \mathrm{d}x \approx \frac{1}{n_f} \sum_{i=1}^{n_f} \frac{f(X_i)g(X_i)w(X_i)}{p_f(X_i)} + \frac{1}{n_g} \sum_{j=1}^{n_g} \frac{f(X_j)g(X_j)w(X_j)}{p_g(X_j)}, \end{equation}

where $n_f$ is the number of samples drawn proportional to $p_f$, $n_g$ the number of samples drawn proportional to $p_g$, and $w_f$ & $w_g$ are special normalization weighting functions chosen so as to not bias the value of the estimator.

Many such weighting functions exist, and a provably good choice is based on the balance heuristic, given by the following equation, where $s \in \langle f,g \rangle = I$ :

\begin{equation} w_s(x) = \frac{n_s\, p_s(x)}{\sum_{i \in I} n_i\, p_i(x)}, \end{equation}

Your task is to implement MIS with this heuristic.

Sampling Emitters and BRDFs (40 points)

Implement DirectMISIntegrator::Li() in src/integrators/directMIS.cpp. This class takes the two measure used for importance sampling measure1 and measure2, the number of emitter samples nSamples1 (in the context of measure1 : solid-angle), the number of BRDF samples nSamples2 (in the context of measure2 : bsdf) and the used weighting functions heuristic as properties which can be retrieved in an XML scene file as usual:

<!-- Use a direct illumination integrator -->
<integrator type="direct-mis">
    <integer name="nSamples1" value="1"/>
    <integer name="nSamples2" value="1"/>
    <string name="measure1" value="solid-angle"/>
    <string name="measure2" value="bsdf"/>
    <string name="heuristic" value="balance"/>
</integrator>

Implementation-wise, your DirectMISIntegrator::Li() function should have two Monte Carlo loops: one over nSamples1 and one over nSamples2. The former should be reminiscent of Assignment 2, except you will need to adapt it to take the weighting function into account. Use the function MIS::balanceHeuristic() in src/warp/mis.cpp for this purpose. For simplicity, you can assume solid angle emitter sampling for this part. The BRDF sampling part isn't much different: once you find an intersection its, you'll want to call BSDF::sample() instead of BSDF::eval() to sample a direction according to the shading point's BRDF.

Rendering the Veach Scene (10 points)

Implement the MIS integrator and render the scene scenes/hw3/Tests/veach-mis.xml containing the popularized test scene from Eric Veach's thesis. Below is a reference image that was rendered to convergence, for comparison.

Testing your Implementation

Implementing MIS can be tricky and you might find many subtle bugs in your previous implementations that didn't surface until now. Here are a few tips to assist you in debugging your code:

• Test with diffuse materials to begin with. Since src/bsdfs/diffuse.cpp is provided in the base code, its Diffuse::sample() method is already correct. Only test MIS with Phong once you're convinced your integrator works properly with diffuse surfaces.
• Test with emitter samples only by setting the number of BSDF samples to zero, and make sure the output image is the same as in Assignment 2.
• Then test with BSDF samples only by setting the number of emitter samples to zero. The output image should obviously have more noise, but should eventually converge with enough samples.
• Finally, combine the two and compare the rendered image to the provided test scene.
• When in absolute doubt, recall that a naive uniform sampling MC estimator will (eventually) generate the correct answer: this is one way to isolate bugs in your specialized PDF sampling and evaluation code, since uniform sampling does not rely on much more than evaluating the integrand terms (as opposed to drawing samples according to them).

---

What to submit

Finished? Render the scenes contained in scenes/hw3/Finals and use the given script to submit your files as you did for the previous assignments. Include any specific comments in the appropriate section of the script.