Subsurface scattering (or SSS) is a mechanism of light transport in which light penetrates the surface of a translucent object, is scattered by interacting with the material, and exits the surface at a different point. Subsurface scattering is important in 3D computer graphics, being necessary for the realistic rendering of materials such as marble, skin and milk. However, SSS is also one of the most complex materials in the state-of-the-art rendering. One reason for this is that the radiance at each point of the object depends on all the other points' radiance and the object's optical thickness, which makes the rendering equation of SSS material a complicated differential equation. So the computational cost of solving this equation is significant, let alone the varying-density property of the material.

Jensen's paper <A Practical Model for Subsurface Light Transport> proposed a robust SSS model, which is a

But instead of evaluating the whole BSSRDF when rendering, what I implemented is still the Single-Scattering Integrator in PBRT, which only counts for self-emission and incident radiance due to direct illumination, incident radiance due to multiple scattering ignored. But the result was nice enough to approximate the whole model.

Let's first take a look at the integral equation of light transfer.

*Bidirectional surface scattering reflectance distribution function (BSSRDF)*. However when in the first contact, those equations in this paper are really beyond the range of my interpretation. But after reading the book <Physically Based Rendering> about Volume Scattering, I got through and understood Jensen's paper.But instead of evaluating the whole BSSRDF when rendering, what I implemented is still the Single-Scattering Integrator in PBRT, which only counts for self-emission and incident radiance due to direct illumination, incident radiance due to multiple scattering ignored. But the result was nice enough to approximate the whole model.

Let's first take a look at the integral equation of light transfer.

At the first glance, this complex equation may be so complicated that it makes no sense at all. But what it really says is simply that the overall differential change in radiance at a point p' along a ray is found by adding the radiance attenuated by σt, which accounts for all the processes that reduce energy along the ray (absorbtion and out-scattering), and the source term (S), which accounts for processes that increase energy along the ray. Given this, we can easily derive the general integral equation that describes light transfer.

Since we only are interested in condition that rays are never blocked inside the material, the first term of this equation is omitted. And in this integral equation, Tr represents the beam transmittance between point p0 and p, which gives the fraction of radiance transmitted between two points on a ray:

where d represents the distance between p and p', w is the normalized direction vector between them. In a homogeneous medium, σt is always a constant, therefore this integral can be trivially evaluated by:

The definition of the source term(S) is given below:

The source term accounts for all the processes that add energy to a ray by summing the emission (Lve) and in-scattering. To get in-scattering, we need to perform a spherical integral at the point p. The phase function (p) is the volumetric analog to BSDF, defining probability distributions for scattering in a particular direction. And the integral is multiplied by a scattering coefficient σs.

Since we need to evaluate the source term (S) at every point along the ray, which require an integral around a sphere, this process is very computationally expensive. However, if we ignore radiance due to multiple scattering, and only take the incident radiance due to direct illumination from each light source into account, the source term can be vastly simplified to:

Since we need to evaluate the source term (S) at every point along the ray, which require an integral around a sphere, this process is very computationally expensive. However, if we ignore radiance due to multiple scattering, and only take the incident radiance due to direct illumination from each light source into account, the source term can be vastly simplified to:

Instead of performing integral around a sphere to get the incidence radiance at point p, we simply sum up all the attenuated radiance from each light source. This procedure can be illustrated by the the figure below:

And since phase functions are normalized, so the isotropic phase function is the constant 1/4π.

With all the simplification above, now the integral equation of light transfer can be trivially written as:

With all the simplification above, now the integral equation of light transfer can be trivially written as:

This is already a huge simplification to the original transfer equation. We can simply use the Monte Carlo method to solve this integral.

In practice, this method is well suited for ray tracing since we evaluated the radiance transfer along a ray. In order to render SSS material, we first evaluate the color of the reflected radiance as normal, and then we trace the ray further into the surface and find the exit point of the ray. Once we get the fragment along which the ray passes through the SSS media, we can evaluate the radiance due to scattering by solving the integral above.

The result can be found

In practice, this method is well suited for ray tracing since we evaluated the radiance transfer along a ray. In order to render SSS material, we first evaluate the color of the reflected radiance as normal, and then we trace the ray further into the surface and find the exit point of the ray. Once we get the fragment along which the ray passes through the SSS media, we can evaluate the radiance due to scattering by solving the integral above.

The result can be found

__here__. I implemented it along with Spherical Harmonics: using the algorithm above to pre-compute the lighting at each vertex. Since it is only convenient for Spherical Harmonics to encode view-independent information, the view-dependent special highlight are evaluated separately in the pixel shader on the fly. So the specular light is not global, which means it will not cast shadows or reflect between surfaces. But the artifact this yields are hardly noticeable and the image produced is already pretty amazing.