## 11.2 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 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 isotropic 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 abstract base class defines the interface for phase functions in pbrt.

<<Media Declarations>>=
class PhaseFunction { public: <<PhaseFunction Interface>>
virtual Float p(const Vector3f &wo, const Vector3f &wi) const = 0; virtual Float Sample_p(const Vector3f &wo, Vector3f *wi, const Point2f &u) const = 0;
};

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.11).

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

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 PhaseHG() function implements this computation.

<<Media Inline Functions>>=
inline Float PhaseHG(Float cosTheta, Float g) { Float denom = 1 + g * g + 2 * g * cosTheta; return Inv4Pi * (1 - g * g) / (denom * std::sqrt(denom)); }

Figure 11.12 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.13 to compare the visual effect of forward- and back-scattering.

HenyeyGreenstein provides a PhaseFunction implementation of the Henyey–Greenstein model.

<<HenyeyGreenstein Declarations>>=
class HenyeyGreenstein : public PhaseFunction { public: <<HenyeyGreenstein Public Methods>>
HenyeyGreenstein(Float g) : g(g) { } Float p(const Vector3f &wo, const Vector3f &wi) const; Float Sample_p(const Vector3f &wo, Vector3f *wi, const Point2f &sample) const;
private: const Float g; };

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

<<HenyeyGreenstein Method Definitions>>=
Float HenyeyGreenstein::p(const Vector3f &wo, const Vector3f &wi) const { return PhaseHG(Dot(wo, wi), g); }

The asymmetry parameter in the Henyey–Greenstein model has a precise meaning. It is the average value of the product of the phase function being approximated and the cosine of the angle between and . 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 aren’t 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 isn’t provided in pbrt but would be easy to add.