Subsurface Scattering

Modified Files:

include/bssrdf.h

include/bsdf.h

src/bssrdf.cpp

src/subsurface.cpp

src/path_mis_sub.cpp

src/path_mis_all.cpp

Implementation

The implementation follows the Sampling Subsurface Reflection Functions and Subsurface Scattering Using the Diffusion Equation of PBRT. The implementation is based on photon beam diffusion technique.
A TabulatedBSSRDF class and a Subsurface class which is a subclass of BSDF is created for the purpose.

BSSRDF Object

The Tabulated BSSRDF is an implementation of a separated version of the bidirectional scattering surface reflectance distribution function, which separates the function into three parts: Fresnel term, profile term, and directional term. The position profile term can be further simplified into a distance profile term under the assumption that only the distance between two points affects the BSSRDF value instead of the specific position. The Fresnel term has a closed form and the directional term is a scaled version of the fresnel transmittance, and the scale is the first moment of the fresnel function with closed formula.
The problem lies in computing and sampling the distance profile term. The Tabulated BSSRDF holds the assumption that the refraction index and phase function parameter stays the same for the whole material. And in my implementation, the other parameters including extinction coefficient and albedo are also fixed for a single material for simplicity. And the photon beam diffusion technique is used here to build the table which holds the pre-computed value of the distance profile term for a range of distance(radius) and albedo values which can be used to evaluate and sample the distance profile term fast during rendering. The PBD method holds several more assumption to compute the scattering term, it assumes the material is homogeneous with highly scattering coefficient to use diffusion approximation which makes the complex equation of light transfer into a relatively simple diffusion equation, and assume the medium is semi-infinite under a plane and place a light source at some depth under the plane to evaluate the radiant exitance at some position on the plane. Based on these, the profile term is evaluated for both multiple scattering and single scattering.
Above we talked about evaluating and sampling the simplified profile function which only considers the distance, but we need to sample an actual point on the surface during rendering. And this is done by sampling direction on a planar approximation of the surface and constructing a ray to intersect with the object. To make it more robust, different axes are sampled to define the surface. All the possible intersections with the object will be found and one is randomly sampled for use.

Subsurface Material

The subsurface material is like any other BSDF but holds a BSSRDF object that shows it requires handling the subsurface scattering. And the BSDF part of the subsurface behaves like a dielectric material.

Integrator Behavior

To handle the subsurface in the path tracer, a branch is added after the emitter sampling and the bsdf sampling. If the BSDF of the intersection object is subsurface BSDF and the sampled ray of the BSDF sampling is a refraction ray, then the subsurface handling is activated.
Firstly, an incident position on the subsurface object is sampled. We also got the pdf of sampling and the scattering factor between the sampled position and the intersection position, which needs to be multiplied to the throughput. Then emitter sampling with MIS is done like before, except the bsdf here is a special one. The BSDF on the sampled position is a diffuse BSDF with the value of the directional term of the BSSRDF. And then, a new direction is sampled with a cosine weighted hemisphere sampling at the sampled position, and the bsdf value is evaluated like in the emitter sampling. Now the ray and the throughput are updated. Remember to set the is_last_specular boolean value to false, as the subsurface scattering is not a specular one.

Validation

For this part, the result is compared with PBRT because Mitsuba doesn't support subsurface scattering yet.

(1) First scene is a simple cbox scene with two spheres with different subsurface materials. All the result is rendered with 256spp.

The scene file are in scenes/subsurface/cbox.

Cbox with two subsurface sphere.

(2) To show the special translucent property of subsurface material, I made a scene with the teapot body and a small but bright sphere emitter inside the teapot to light up the teapot. And we can see the beautiful translucent appearance of the teapot. In the comparison, my result with 256spp has much more noise than PBRT, and it becomes better with 512spp. And I found that my rendering with subsurface scattering is much faster than PBRT with the same spp: Mine with 256spp takes ~16s, with 512spp takes ~32s and PBRT with 256spp takes ~26s.

The scene file are in scenes/subsurface/teapot.

Teapot with a emitter inside.

And due to the noise, it's not very obvious that my rendering has the same distribution with PBRT's. Therefore, I show a comparison of the two in Tev. As the screenshot shows, the mean error is 0.

Comparison in Tev