## 10.1 Sampling and Antialiasing

The sampling task from Chapter 7 was a
frustrating one since the aliasing problem was known to be unsolvable from
the start. The infinite frequency content of geometric edges and hard
shadows *guarantees* aliasing in the final images, no matter how high
the image sampling rate. (Our only consolation is that the visual impact
of this remaining aliasing can be reduced to unobjectionable levels with a
sufficient number of well-placed samples.)

Fortunately, for textures things are not this difficult from the start: either there is often a convenient analytic form of the texture function available, which makes it possible to remove excessively high frequencies before sampling it, or it is possible to be careful when evaluating the function so as not to introduce high frequencies in the first place. When this problem is carefully addressed in texture implementations, as is done through the rest of this chapter, there is usually no need for more than one sample per pixel in order to render an image without texture aliasing.

Two problems must be addressed in order to remove aliasing from texture functions:

- The sampling rate in texture space must be computed. The screen space sampling rate is known from the image resolution and pixel sampling rate, but here we need to determine the resulting sampling rate on a surface in the scene in order to find the rate at which the texture function is being sampled.
- Given the texture sampling rate, sampling theory must be applied to guide the computation of a texture value that doesn’t have higher frequency variation than can be represented by the sampling rate (e.g., by removing excess frequencies beyond the Nyquist limit from the texture function).

These two issues will be addressed in turn throughout the rest of this section.

### 10.1.1 Finding the Texture Sampling Rate

Consider an arbitrary texture function that is a function of position, , defined on a surface in the scene. If we ignore the complications introduced by visibility issues—the possibility that another object may occlude the surface at nearby image samples or that the surface may have a limited extent on the image plane—this texture function can also be expressed as a function over points on the image plane, , where is the function that maps image points to points on the surface. Thus, gives the value of the texture function as seen at image position .

As a simple example of this idea, consider a 2D texture function applied to a quadrilateral that is perpendicular to the axis and has corners at the world space points , , , and . If an orthographic camera is placed looking down the axis such that the quadrilateral precisely fills the image plane and if points on the quadrilateral are mapped to 2D texture coordinates by

then the relationship between and screen pixels is straightforward:

where the overall image resolution is (Figure 10.2). Thus, given a sample spacing of one pixel in the image plane, the sample spacing in texture parameter space is , and the texture function must remove any detail at a higher frequency than can be represented at that sampling rate.

This relationship between pixel coordinates and texture coordinates, and thus the relationship between their sampling rates, is the key bit of information that determines the maximum frequency content allowable in the texture function. As a slightly more complex example, given a triangle with texture coordinates at the vertices and viewed with a perspective projection, it’s possible to analytically find the differences in and across the sample points on the image plane. This is the basis of basic texture map antialiasing in graphics processors.

For more complex scene geometry, camera projections, and mappings to texture coordinates, it is much more difficult to precisely determine the relationship between image positions and texture parameter values. Fortunately, for texture antialiasing, we don’t need to be able to evaluate for arbitrary but just need to find the relationship between changes in pixel sample position and the resulting change in texture sample position at a particular point on the image. This relationship is given by the partial derivatives of this function, and . For example, these can be used to find a first-order approximation to the value of ,

If these partial derivatives are changing slowly with respect to the distances and , this is a reasonable approximation. More importantly, the values of these partial derivatives give an approximation to the change in texture sample position for a shift of one pixel in the and directions, respectively, and thus directly yield the texture sampling rate. For example, in the previous quadrilateral example, , , , and .

The key to finding the values of these partial derivatives in the general
case lies in the `RayDifferential` structure, which was defined in
Section 2.5.1. This structure is initialized for
each camera ray by the `Camera::GenerateRayDifferential()` method; it
contains not only the ray actually being traced through the scene but also
two additional rays, one offset horizontally one pixel sample from the camera ray
and the other offset vertically by one pixel sample. All of the geometric ray
intersection routines use only the main camera ray for their computations;
the auxiliary rays are ignored (this is easy to do because
`RayDifferential` is a subclass of `Ray`).

Here we will use the offset rays to estimate the partial
derivatives of the mapping p from image position to world space
position and the partial derivatives of the mappings and from to
parametric coordinates, giving the partial derivatives of world space positions
and and the partial derivatives of parametric coordinates
, , , and . In Section 10.2, we will see
how these can be used to compute the screen space
derivatives of arbitrary quantities based on or and
consequently the sampling rates of these quantities. The values of these
partial derivatives at the intersection point are stored in the
`SurfaceInteraction` structure. They are declared as
`mutable`, since they are set in a method that takes a `const`
instance of that object.

The `SurfaceInteraction::ComputeDifferentials()` method computes these
values. It is called by
`SurfaceInteraction::ComputeScatteringFunctions()` before the
`Material`’s `ComputeScatteringFunctions()` method is called so that
these values will be available for any texture evaluation routines that are
called by the material. Because ray differentials aren’t available for all
rays traced by the system (e.g., rays starting from light sources traced
for photon mapping or bidirectional path tracing), the
`hasDifferentials` field of the `RayDifferential` must be checked
before these computations are performed. If the differentials are not present,
then the derivatives are all set to zero (which will eventually lead
to unfiltered point sampling of textures).

`A`,

`Bx`, and

`By`matrices for offset computation>>

The key to computing these estimates is the assumption that the surface is locally flat with respect to the sampling rate at the point being shaded. This is a reasonable approximation in practice, and it is hard to do much better. Because ray tracing is a point-sampling technique, we have no additional information about the scene in between the rays we have traced. For highly curved surfaces or at silhouette edges, this approximation can break down, though this is rarely a source of noticeable error in practice.

For this approximation, we need the plane through the point intersected by the main ray that is tangent to the surface. This plane is given by the implicit equation

where , , , and . We can then compute the intersection points and between the auxiliary rays and and this plane (Figure 10.3). These new points give an approximation to the partial derivatives of position on the surface and , based on forward differences:

Because the differential rays are offset one pixel sample in each direction, there’s no need to divide these differences by a value, since .

`A`,

`Bx`, and

`By`matrices for offset computation>>

The ray–plane intersection algorithm described in Section 3.1.2 gives the value where a ray described by origin and direction intersects a plane described by :

To compute this value for the two auxiliary rays, the
plane’s coefficient is computed first. It isn’t necessary to compute the , , and
coefficients, since they’re available in `n`. We can then apply
the formula directly.

Using the positions and , an approximation to their respective coordinates can be found by taking advantage of the fact that the surface’s partial derivatives and form a (not necessarily orthogonal) coordinate system on the plane and that the coordinates of the auxiliary intersection points in terms of this coordinate system are their coordinates with respect to the parameterization (Figure 10.4).

Given a position on the plane, we can compute its position with respect to the coordinate system by

or, equivalently,

A solution to this linear system of equations for one of the auxiliary points or gives the corresponding screen space partial derivatives or , respectively.

This linear system has three equations with two unknowns—that is, it’s
overconstrained. We need to be careful since one of the equations may be
degenerate—for example, if and are in the plane such
that their components are both zero, then the third equation will be
degenerate. Therefore, we’d like to solve the system of equations using
two equations that don’t give a degenerate system. An easy way to do
this is to take the cross product of and , see which
coordinate of the result has the largest magnitude, and use the other two.
Their cross product is already available in `n`, so using this
approach is straightforward. Even after all this, it may happen that the
linear system has no solution (usually due to the partial derivatives not
forming a coordinate system on the plane). In that case, all that can be
done is to return arbitrary values.

`A`,

`Bx`, and

`By`matrices for offset computation>>

`A`,

`Bx`, and

`By`matrices for offset computation>>=

### 10.1.2 Filtering Texture Functions

It is necessary to
remove frequencies in texture functions that are past the Nyquist limit for
the texture sampling rate. The goal is to compute, with as few approximations as
possible, the result of the *ideal texture resampling* process, which says
that in order to evaluate without aliasing, we must first
band-limit it, removing frequencies beyond the Nyquist limit by
convolving it with the sinc filter:

The band-limited function in turn should then be convolved with the pixel filter centered at the point on the screen at which we want to evaluate the texture function:

This gives the theoretically perfect value for the texture as projected onto the screen.

In practice, there are many simplifications that can be made to this process, with little reduction in visual quality. For example, a box filter may be used for the band-limiting step, and the second step is usually ignored completely, effectively acting as if the pixel filter were a box filter, which makes it possible to do the antialiasing work completely in texture space and simplifies the implementation significantly. The EWA filtering algorithm in Section 10.4.5 is a notable exception in that it assumes a Gaussian pixel filter.

The box filter is easy to use, since it can be applied
analytically by computing the average of the texture function over the
appropriate region. Intuitively, this is a reasonable approach to the
texture filtering problem, and it can be computed directly for many texture
functions. Indeed, through the rest of this chapter, we will often use a
box filter to average texture function values between samples and
informally use the term *filter region* to describe the area being
averaged over. This is the most common approach when filtering texture
functions.

Even the box filter, with all of its shortcomings, gives acceptable results for texture filtering in many cases. One factor that helps is the fact that a number of samples are usually taken in each pixel. Thus, even if the filtered texture values used in each one are sub-optimal, once they are filtered by the pixel reconstruction filter, the end result generally doesn’t suffer too much.

An alternative to using the box filter to filter texture functions is to
use the observation that the effect of the ideal sinc filter is to let
frequency components below the Nyquist limit pass through unchanged but to
remove frequencies past it. Therefore, if we know the frequency content of
the texture function (e.g., if it is a sum of terms, each one with known
frequency content), then if we replace the high-frequency terms with their
average values, we are effectively doing the work of the sinc prefilter.
This approach is known as *clamping* and is the basis for antialiasing
in the textures based on the noise function in Section 10.6.

Finally, for texture functions where none of these techniques is easily
applied, a final option is *supersampling*—the function
is evaluated and filtered at multiple locations near the main evaluation
point, thus increasing the sampling rate in texture space. If a box filter
is used to filter these sample values, this is equivalent to averaging the
value of the function. This approach can be expensive if the texture
function is complex to evaluate, and as with image sampling a very large
number of samples may be needed to remove aliasing. Although this is a
brute-force solution, it is still more efficient than increasing the image
sampling rate, since it doesn’t incur the cost of tracing more rays through
the scene.

### 10.1.3 Ray Differentials for Specular Reflection and Transmission

Given the effectiveness of ray differentials for finding filter regions for
texture antialiasing for camera rays, it is useful to extend the method
to make it possible to determine texture space sampling rates for objects
that are seen indirectly via specular reflection or refraction; objects
seen in mirrors, for example, should also no more have texture aliasing
than directly visible objects. Igehy (1999) developed an elegant solution to
the problem of how to find the appropriate differential rays for specular
reflection and refraction, which is the approach used in
`pbrt`.

Figure 10.5 illustrates the difference that proper texture filtering for specular reflection and transmission can make: it shows a glass ball and a mirrored ball on a plane with a texture map containing high-frequency components. Ray differentials ensure that the images of the texture seen via reflection and refraction from the balls are free of aliasing artifacts. Here, ray differentials eliminate aliasing without excessively blurring the texture.

In order to compute the reflected or transmitted ray differentials at a surface intersection point, we need an approximation to the rays that would have been traced at the intersection points for the two offset rays in the ray differential that hit the surface (Figure 10.6). The new ray for the main ray is computed by the BSDF, so here we only need to compute the outgoing rays for the and differentials.

`pbrt`estimates the change in reflected direction as a function of image space position and approximates the ray differential’s direction with the main ray’s direction added to the estimated change in direction.

For both reflection and refraction, the origin of each differential ray is
easily found. The `SurfaceInteraction::ComputeDifferentials()`
method previously computed approximations for how much the surface position
changes with respect to position on the image plane and
. Adding these offsets to the intersection point of the main ray
gives approximate origins for the new rays. If the incident ray doesn’t
have differentials, then it’s impossible to compute reflected ray
differentials and this step is skipped.

`rd`for specular reflection>>=

Finding the directions of these rays is slightly trickier. Igehy (1999) observed that if we know how much the reflected direction changes with respect to a shift of a pixel sample in the and directions on the image plane, we can use this information to approximate the direction of the offset rays. For example, the direction for the ray offset in is

Recall from Equation (8.5) that for a general world space surface normal and outgoing direction, the direction for perfect specular reflection is

Fortunately, the partial derivatives of this expression are easily computed:

Using the properties of the dot product, it can be shown that

The value of can be found from the difference
between the direction of the ray differential’s main ray and the direction
of the offset ray, and
all of the other necessary quantities are readily available from the
`SurfaceInteraction`, so the implementation of this computation for
the partial derivatives in and is straightforward.

A similar process of differentiating the equation for the direction of a
specularly transmitted ray, Equation (8.8), gives
the equation to find the differential change in the transmitted direction.
We won’t include the derivation or our implementation here, but refer the
interested reader to the original paper and to the `pbrt` source code,
respectively.