11.3 Phase Functions

Just as there is a wide variety of BSDF models that describe scattering from surfaces, many phase functions have also been developed. These range from parameterized models (which can be used to fit a function with a small number of parameters to measured data) to analytic models that are based on deriving the scattered radiance distribution that results from particles with known shape and material (e.g., spherical water droplets).

In most naturally occurring media, the phase function is a 1D function of the angle theta between the two directions omega Subscript normal o and omega Subscript normal i ; these phase functions are often written as p left-parenthesis cosine theta right-parenthesis . Media with this type of phase function are called isotropic or symmetric because their response to incident illumination is (locally) invariant under rotations. In addition to being normalized, an important property of naturally occurring phase functions is that they are reciprocal: the two directions can be interchanged and the phase function’s value remains unchanged. Note that symmetric phase functions are trivially reciprocal because cosine left-parenthesis negative theta right-parenthesis equals cosine left-parenthesis theta right-parenthesis .

In anisotropic media that consist of particles arranged in a coherent structure, the phase function can be a 4D function of the two directions, which satisfies a more involved kind of reciprocity relation. Examples of this are crystals or media made of coherently oriented fibers; the “Further Reading” discusses these types of media further.

In a slightly confusing overloading of terminology, phase functions themselves can be isotropic or anisotropic as well. Thus, we might have an anisotropic phase function in an isotropic medium. An isotropic phase function describes equal scattering in all directions and is thus independent of either of the two directions. Because phase functions are normalized, there is only one such function:

p left-parenthesis omega Subscript normal o Baseline comma omega Subscript normal i Baseline right-parenthesis equals StartFraction 1 Over 4 pi EndFraction period

The PhaseFunction class defines the PhaseFunction interface. Only a single phase function is currently provided in pbrt, but we have used the TaggedPointer machinery to make it easy to add others. Its implementation is in the file base/medium.h.

<<PhaseFunction Definition>>= 
class PhaseFunction : public TaggedPointer<HGPhaseFunction> { public: <<PhaseFunction Interface>> 
using TaggedPointer::TaggedPointer; std::string ToString() const; Float p(Vector3f wo, Vector3f wi) const; pstd::optional<PhaseFunctionSample> Sample_p(Vector3f wo, Point2f u) const; Float PDF(Vector3f wo, Vector3f wi) const;
};

The p() method returns the value of the phase function for the given pair of directions. As with BSDFs, pbrt uses the convention that the two directions both point away from the point where scattering occurs; this is a different convention from what is usually used in the scattering literature (Figure 11.13).

<<PhaseFunction Interface>>= 
Float p(Vector3f wo, Vector3f wi) const;

Figure 11.13: Phase functions in pbrt are implemented with the convention that both the incident direction and the outgoing direction point away from the point where scattering happens. This is the same convention that is used for BSDFs in pbrt but is different from the convention in the scattering literature, where the incident direction generally points toward the scattering point. The angle between the two directions is denoted by  theta .

It is also useful to be able to draw samples from the distribution described by a phase function. PhaseFunction implementations therefore must provide a Sample_p() method, which samples an incident direction omega Subscript normal i given the outgoing direction omega Subscript normal o and a sample value in left-bracket 0 comma 1 right-parenthesis squared .

<<PhaseFunction Interface>>+=  
pstd::optional<PhaseFunctionSample> Sample_p(Vector3f wo, Point2f u) const;

Phase function samples are returned in a structure that stores the phase function’s value p, the sampled direction wi, and the PDF pdf.

<<PhaseFunctionSample Definition>>= 
struct PhaseFunctionSample { Float p; Vector3f wi; Float pdf; };

An accompanying PDF() method returns the value of the phase function sampling PDF for the provided directions.

<<PhaseFunction Interface>>+= 
Float PDF(Vector3f wo, Vector3f wi) const;

11.3.1 The Henyey–Greenstein Phase Function

A widely used phase function was developed by Henyey and Greenstein (1941). This phase function was specifically designed to be easy to fit to measured scattering data. A single parameter g (called the asymmetry parameter) controls the distribution of scattered light:

p Subscript normal upper H normal upper G Baseline left-parenthesis cosine theta right-parenthesis equals StartFraction 1 Over 4 pi EndFraction StartFraction 1 minus g squared Over left-parenthesis 1 plus g squared plus 2 g left-parenthesis cosine theta right-parenthesis right-parenthesis Superscript 3 slash 2 Baseline EndFraction period

The HenyeyGreenstein() function implements this computation.

<<Scattering Inline Functions>>+= 
Float HenyeyGreenstein(Float cosTheta, Float g) { Float denom = 1 + Sqr(g) + 2 * g * cosTheta; return Inv4Pi * (1 - Sqr(g)) / (denom * SafeSqrt(denom)); }

The asymmetry parameter g in the Henyey–Greenstein model has a precise meaning. It is the integral of the product of the given phase function and the cosine of the angle between omega prime Subscript and omega Subscript and is referred to as the mean cosine. Given an arbitrary phase function p , the value of g can be computed as

g equals integral Underscript script upper S squared Endscripts p left-parenthesis minus omega Subscript Baseline comma omega prime Subscript Baseline right-parenthesis left-parenthesis omega Subscript Baseline dot omega prime Subscript Baseline right-parenthesis normal d omega Subscript Baseline Superscript prime Baseline equals 2 pi integral Subscript 0 Superscript pi Baseline p left-parenthesis minus cosine theta right-parenthesis cosine theta sine theta normal d theta Subscript Baseline period

Thus, an isotropic phase function gives g equals 0 , as expected.

Any number of phase functions can satisfy this equation; the g value alone is not enough to uniquely describe a scattering distribution. Nevertheless, the convenience of being able to easily convert a complex scattering distribution into a simple parameterized model is often more important than this potential loss in accuracy.

More complex phase functions that are not described well with a single asymmetry parameter can often be modeled by a weighted sum of phase functions like Henyey–Greenstein, each with different parameter values:

p left-parenthesis omega Subscript Baseline comma omega prime Subscript Baseline right-parenthesis equals sigma-summation Underscript i equals 1 Overscript n Endscripts w Subscript i Baseline p Subscript i Baseline left-parenthesis omega Subscript Baseline right-arrow omega prime Subscript Baseline right-parenthesis comma

where the weights w Subscript i sum to one to maintain normalization. This generalization is not provided in pbrt but would be easy to add.

Figure 11.14 shows plots of the Henyey–Greenstein phase function with varying asymmetry parameters. The value of g for this model must be in the range left-parenthesis negative 1 comma 1 right-parenthesis . Negative values of g correspond to back-scattering, where light is mostly scattered back toward the incident direction, and positive values correspond to forward-scattering. The greater the magnitude of g , the more scattering occurs close to the omega Subscript or minus omega Subscript directions (for back-scattering and forward-scattering, respectively). See Figure 11.15 to compare the visual effect of forward- and back-scattering.

Figure 11.14: Plots of the Henyey–Greenstein Phase Function for Asymmetry g Parameters minus 0.25 and 0.7. Negative g values describe phase functions that primarily scatter light back in the incident direction, and positive g values describe phase functions that primarily scatter light forward in the direction it was already traveling (here, along the plus x axis).

Figure 11.15: Ganesha model filled with participating media rendered with (left) strong backward scattering ( g equals negative 0.9 ) and (right) strong forward scattering ( g equals 0.9 ). Because most of the light comes from a light source behind the objects, forward scattering leads to more light reaching the camera in this case.

The HGPhaseFunction class implements the Henyey–Greenstein model in the context of the PhaseFunction interface.

<<HGPhaseFunction Definition>>= 
class HGPhaseFunction { public: <<HGPhaseFunction Public Methods>> 
HGPhaseFunction(Float g) : g(g) {} Float p(Vector3f wo, Vector3f wi) const { return HenyeyGreenstein(Dot(wo, wi), g); } pstd::optional<PhaseFunctionSample> Sample_p(Vector3f wo, Point2f u) const { Float pdf; Vector3f wi = SampleHenyeyGreenstein(wo, g, u, &pdf); return PhaseFunctionSample{pdf, wi, pdf}; } Float PDF(Vector3f wo, Vector3f wi) const { return p(wo, wi); } static const char *Name() { return "Henyey-Greenstein"; } std::string ToString() const;
private: <<HGPhaseFunction Private Members>> 
Float g;
};

Its only parameter is g , which is provided to the constructor and stored in a member variable.

<<HGPhaseFunction Public Methods>>= 
HGPhaseFunction(Float g) : g(g) {}

<<HGPhaseFunction Private Members>>= 
Float g;

Evaluating the phase function is a simple matter of calling the HenyeyGreenstein() function.

<<HGPhaseFunction Public Methods>>+=  
Float p(Vector3f wo, Vector3f wi) const { return HenyeyGreenstein(Dot(wo, wi), g); }

It is possible to sample directly from the Henyey–Greenstein phase function’s distribution. This operation is provided via a stand-alone utility function. Because the sampling algorithm is exact and because the Henyey–Greenstein phase function is normalized, the PDF is equal to the phase function’s value for the sampled direction.

<<Sampling Function Definitions>>+=  
Vector3f SampleHenyeyGreenstein(Vector3f wo, Float g, Point2f u, Float *pdf) { <<Compute cosine theta for Henyey–Greenstein sample>> 
Float cosTheta; if (std::abs(g) < 1e-3f) cosTheta = 1 - 2 * u[0]; else cosTheta = -1 / (2 * g) * (1 + Sqr(g) - Sqr((1 - Sqr(g)) / (1 + g - 2 * g * u[0])));
<<Compute direction wi for Henyey–Greenstein sample>> 
Float sinTheta = SafeSqrt(1 - Sqr(cosTheta)); Float phi = 2 * Pi * u[1]; Frame wFrame = Frame::FromZ(wo); Vector3f wi = wFrame.FromLocal(SphericalDirection(sinTheta, cosTheta, phi));
if (pdf) *pdf = HenyeyGreenstein(cosTheta, g); return wi; }

The PDF for the Henyey–Greenstein phase function is separable into theta and phi components, with p left-parenthesis phi right-parenthesis equals 1 slash left-parenthesis 2 pi right-parenthesis as usual. The main task is to sample cosine theta . With pbrt’s convention for the orientation of direction vectors, the distribution for theta is

cosine theta equals minus StartFraction 1 Over 2 g EndFraction left-parenthesis 1 plus g squared minus left-parenthesis StartFraction 1 minus g squared Over 1 plus g minus 2 g xi Subscript Baseline EndFraction right-parenthesis squared right-parenthesis

if g not-equals 0 ; otherwise, cosine theta equals 1 minus 2 xi Subscript gives a uniform sampling over the sphere of directions.

<<Compute cosine theta for Henyey–Greenstein sample>>= 
Float cosTheta; if (std::abs(g) < 1e-3f) cosTheta = 1 - 2 * u[0]; else cosTheta = -1 / (2 * g) * (1 + Sqr(g) - Sqr((1 - Sqr(g)) / (1 + g - 2 * g * u[0])));

The left-parenthesis cosine theta comma phi right-parenthesis values specify a direction with respect to a coordinate system where wo is along the plus z axis. Therefore, it is necessary to transform the sampled vector to wo’s coordinate system before returning it.

<<Compute direction wi for Henyey–Greenstein sample>>= 
Float sinTheta = SafeSqrt(1 - Sqr(cosTheta)); Float phi = 2 * Pi * u[1]; Frame wFrame = Frame::FromZ(wo); Vector3f wi = wFrame.FromLocal(SphericalDirection(sinTheta, cosTheta, phi));

The HGPhaseFunction sampling method is now easily implemented.

<<HGPhaseFunction Public Methods>>+=  
pstd::optional<PhaseFunctionSample> Sample_p(Vector3f wo, Point2f u) const { Float pdf; Vector3f wi = SampleHenyeyGreenstein(wo, g, u, &pdf); return PhaseFunctionSample{pdf, wi, pdf}; }

Because sampling is exact and phase functions are normalized, its PDF() method just evaluates the phase function for the given directions.

<<HGPhaseFunction Public Methods>>+= 
Float PDF(Vector3f wo, Vector3f wi) const { return p(wo, wi); }