## 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 between the two directions and ; these phase
functions are often written as . 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 .

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:

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`.

The `p()` method returns the value of the phase function for the given
pair of directions. As with `BSDF`s, `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).

`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 .

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 given the outgoing direction and a sample
value in .

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

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

### 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 (called the *asymmetry
parameter*)
controls the distribution of scattered light:

The `HenyeyGreenstein()` function implements this computation.

The asymmetry parameter 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 and and is referred to
as the *mean cosine*. Given an
arbitrary phase function , the value of can be computed as

Thus, an isotropic phase function gives , as expected.

Any number of phase functions can satisfy this equation; the 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:

where the weights 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 for this
model must be in the range . Negative values of 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 , the more scattering occurs close to the
or directions (for back-scattering and forward-scattering,
respectively). See Figure 11.15 to compare the visual
effect of forward- and back-scattering.

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

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

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

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.

`wi`for Henyey–Greenstein sample>>

The PDF for the Henyey–Greenstein phase function is separable into
and components, with as usual. The main
task is to sample . With `pbrt`’s convention for the
orientation of direction vectors, the distribution for is

if ; otherwise, gives a uniform sampling over the sphere of directions.

The values specify a direction with respect to a
coordinate system where `wo` is along the axis. Therefore, it is
necessary to transform the sampled vector to `wo`’s coordinate
system before returning it.

`wi`for Henyey–Greenstein sample>>=

The `HGPhaseFunction` sampling method is now easily implemented.

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