### From BlenderWiki

# Subsurface Scattering

Implementation of raytraced subsurface scattering in Cycles.

### Sampling Nearby Points

The BSSRDF [1] assumes the surface is an infinite plane, and sampling nearby points on such a surface would be simple, but in practice we apply it to arbitrary geometry about which we can make few assumptions.

In [2] the proposed solution is to distribute points on the surface before rendering. The problem is that this method does not fit well with our rendering architecture, it means we must used a fixed number of samples, it gives a preprocessing pass and does not work with progressive rendering or more advanced algorithms like adaptive sampling or bidirectional ray tracing. We want to be able to start rendering immediately when lights or materials change. Another problem is the memory usage and rendering of many objects, estimating how many points need to be sampled on them is hard to do beforehand.

For Blender Internal we generated two samples per pixel (front & back), which works well to keep memory usage predictable, but does not give a good distribution for global illumination where the surface may be hit from any direction.

A different method is to sample points in the tangent plane, and project them downwards onto the surface using raytracing [3]. This doesn't work so well in practice however, as a distance from the shading point in the tangent plane would not have the same distance as the projected point. Further note that on e.g. the corners of a cube, this will not be able to find points on an adjacent side of the cube. While this may be considered a corner case, even if the side is slightly slanted, it will receive very few samples.

Instead we use the method described in Appendix B of [4], which uses sampling of random lines between two points on a sphere. The size of the sphere is varied based on the BSSRDF function, which gives us a method to sample points on nearby surfaces within a certain distance. This gives results that are close to [2], which is a method known to work well for film production. Another advantage compared to tangent plane sampling is that this method is symmetric and so works with bidirectional ray tracing.

However this method still has issues. For a single line sample, we can hit zero, one or multiple surfaces. In case we hit multiple surfaces, we pick one of them at random. The big problem is hitting zero surfaces, which happens frequently and gives noisy results. In the case of an infinite plane, half the lines would miss the surface. Besides the noise, another problem is keeping the light contribution normalized, for an infinite plane we can simply multiply by 2 to make up for the missed samples, but for arbitrary geometry this is not possible.

As a solution we keep taking more line samples until we hit a surface, or a limit is reached. If the limit is reached, we simply stay at the current shading point and continue from there. A higher limit will give more accurate result but wasted time on shading points that have few nearby surfaces. We set the limit to 8, which means for an infinite plane 1/256 of the light contribution will stay at the shading point.

### Importance Sampling the BSSRDF

For line sampling we need to generate distance samples according to the BSSRDF. For importance sampling this means integrating and inverting the falloff function, which seems to be impractical to do analytically (maybe there is a simple formulation, but I couldn't find it). So we we'll do this numerically based on lookup tables, which will likely also be faster.

Our BSSRDF uses the reparametrization presented in [2], which has 3 parameters:

- Diffuse mean free path (average scattering distance)
- Index of refraction
- Diffuse Reflectance.

We write this as:

Rd(r, mfp, ior, refl)

Building a lookup table for all possible combinations would be impractical, but we do want to be able to texture these parameters and so we can't build a single lookup table for fixed parameters. So we make some simplications to make this work.

For the mean free path, we can see that the Rd function assumes no particular units, the index of refraction and reflectance are unitless, and there are no other builtin parameters. This means we must be able to factor out the mean free path parameter if the BSSRDF is to be scale independent, and we can indeed do this:

Rd(r, mfp, ior, refl) = Rd(r/mfp, 1, ior, refl)/(mfp^2)

Further we make the simplifaction that the index of refraction is always 1.3, which is the case for typical surfaces rendered with subsurface scattering. This leaves us with just the distance r and reflectance parameters, which makes building a lookup table practical and fast. Making the index of refraction configurable should be possible too but we leave this to another day.

Building a lookup table for importance sampling is a relatively standard operation, numerically integrating and inverting Rd. Because of the line sampling method we use, we are actually sampling multiple overlapping spheres and we need to take this into account.

### References

- [1] A Practical Model for Subsurface Light Transport
- [2] A Rapid Hierarchical Rendering Technique for Translucent Materials
- [3] Implementing a skin BSSRDF (or several...)
- [4] Bidirectional Lightcuts