## 8.8 Image Reconstruction

As discussed in Section 5.4.3, each pixel in the `Film`
computes an estimate of the integral of the product of a filter function
with samples taken from the image function. In
Section 8.1, we saw that sampling theory provides a
mathematical foundation for how this filtering operation should be
performed in order to achieve an antialiased result. We should, in
principle:

- Reconstruct a continuous image function from the set of image samples.
- Prefilter that function to remove any frequencies past the Nyquist limit for the pixel spacing.
- Sample the prefiltered function at the pixel locations to compute the final pixel values.

Because we know that we will be resampling the function at only the pixel locations, it is not necessary to construct an explicit representation of the function. Instead, we can combine the first two steps using a single filter function.

Recall that if the original function had been uniformly sampled at a frequency greater than the Nyquist frequency and reconstructed with the sinc filter, then the reconstructed function in the first step would match the original image function perfectly—quite a feat since we only have point samples. But because the image function almost always will have higher frequencies than could be accounted for by the sampling rate (due to edges, etc.), we chose to sample it nonuniformly, trading off noise for aliasing.

The theory behind ideal reconstruction depends on the samples being
uniformly spaced. While a number of attempts have been made to extend the
theory to nonuniform sampling, there is not yet an accepted approach to
this problem. Furthermore, because the sampling rate is known to be
insufficient to capture the function, perfect reconstruction is not
possible. Recent research in the field of sampling theory has revisited
the issue of reconstruction with the explicit acknowledgment that perfect
reconstruction is not generally attainable in practice. This slight shift
in perspective has led to powerful new reconstruction techniques. In
particular, the goal of research in reconstruction theory has shifted from
perfect reconstruction to developing reconstruction techniques that can be
shown to minimize error between the reconstructed function and the original
function, *regardless of whether the original was band limited*.

The sinc filter is not an appropriate choice here: recall that the ideal
sinc filter is prone to ringing when the underlying function has
frequencies beyond the Nyquist limit, meaning edges in
the image have faint replicated copies of the edge in nearby pixels (the
Gibbs phenomenon; see Section 8.1.5).
Furthermore, the sinc filter has *infinite support*: it does not fall
off to zero at a finite distance from its center, so all the image
samples would need to be filtered for each output pixel. In practice,
there is no single best filter function. Choosing the best one for a
particular scene takes a mixture of quantitative evaluation and qualitative
judgment. `pbrt` therefore provides a variety of choices.

Figure 8.48 shows comparisons of zoomed-in regions of images rendered using a variety of the filters from this section to reconstruct pixel values.

*(Crown model courtesy of Martin Lubich.)*

### 8.8.1 Filter Interface

The `Filter` class defines the interface for pixel reconstruction
filters in `pbrt`. It is defined in the file `base/filter.h`.

All filters are 2D functions that are centered at the origin and define a radius beyond
which they have a value of 0. The radii are different in the and
directions but are assumed to be symmetric in each. A filter provides
its radius via the `Radius()` method. The filter’s overall extent
in each direction (its *support*) is twice the value of its
corresponding radius (Figure 8.49).

`Filter` implementations must also provide a method that evaluates
their filter function. This function may be called with points that are
outside of the filter’s radius; in this case it is the
responsibility of the implementation to detect this case and return the
value 0. It is not required for the filter values returned by
`Evaluate()` to integrate to 1 since the estimator used to compute pixel values,
Equation (5.13), is self-normalizing.

Filters also must be able to return their integral. Most are able to
compute this value in closed form. Thus, if calling code requires a
normalized filter function, it is easy enough to find it by dividing values
returned by `Evaluate()` by the integral.

Filters must also provide an importance sampling method, `Sample`,
which takes a random sample `u` in .

The returned `FilterSample` structure stores both the sampled position
`p` and a weight, which is the ratio of the value of the filter
function at `p` to the value of the PDF used for sampling there.
Because some filters are able to exactly sample from their distribution,
returning this ratio directly allows them to save the trouble of evaluating
those two values and instead to always return a weight of 1.

Given the specification of this interface, we can now implement
the `GetCameraSample()` function that most integrators
use to compute the `CameraSample`s that they pass to the
`Camera::GenerateRay()` methods.

`CameraSample`member variables>> return cs; }

One subtlety in the definition of this function is that it is templated
based on the type of the sampler passed to it. If a value of type
`Sampler` is passed to this method, then it proceeds using `pbrt`’s
usual dynamic dispatch mechanism to call the corresponding methods of the
`Sampler` implementation. However, if a concrete sampler type
(e.g., `HaltonSampler`) is passed to it, then the corresponding methods
can be called directly (and are generally expanded inline in the function).
This capability is used to improve performance in `pbrt`’s GPU rendering
path; see Section 15.3.3.

After the filter’s `Sample()` method has returned a
`FilterSample`, the image sample position can be found by adding the
filter’s sampled offset to the pixel coordinate before a shift of in
each dimension accounts for the mapping from discrete to continuous pixel
coordinates (recall Section 8.1.4). The filter’s weight
is passed along in the `CameraSample` so that it is available to the
`Film` when its `AddSample()` method
is called.

`CameraSample`member variables>>=

### 8.8.2 FilterSampler

Not all `Filter`s are able to easily sample from the distributions of
their filter functions. Therefore, `pbrt` provides a `FilterSampler`
class that wraps up the details of sampling based on a tabularized
representation of the filter.

Only the `Filter` and an allocator are provided to the constructor.
We have not found it particularly useful to allow the caller to specify the
rate at which the filter function is sampled to construct the table used
for sampling, so instead hardcode a sampling rate of 32 times per unit
filter extent in each dimension.

`f`>>

`domain` gives the bounds of the filter and `f` stores
tabularized filter function values.

All the filters currently implemented in `pbrt` are symmetric about the
origin, which means that they could be tabularized over a single
quadrant. Further, they are all separable into the product of two 1D
functions. Either of these properties could be exploited to
reduce the amount of storage required for the tables used for sampling.
However, to allow full flexibility with the definition of additional filter
functions, the `FilterSampler` simply evaluates the filter at
equally spaced positions over its entire domain to initialize the `f`
array.

`f`>>=

Given a tabularized function, it is easy to initialize the sampling distribution.

There are two important details in the implementation of its
`Sample()` method. First, the implementation does not
use `Filter::Evaluate()` to evaluate the filter function but instead
uses the tabularized version of it in `f`. By using the piecewise
constant approximation of the filter function, it ensures that the returned
weight for a sampled point is always for a
constant . If it did not do this, there
would be variation in the returned weight for non-constant filter functions,
due to the sampling distribution not being exactly proportional to the
filter function—see Figure 8.50, which
illustrates the issue.

`FilterSampler`. If filter positions are found by sampling from and contributions are weighted using the ratio , then different samples may have very different contributions. For example, the two points shown have a difference in their values. This variation in filter weights can lead to variance in images and therefore the

`FilterSampler`uses the same piecewise-constant approximation of for evaluation as is used for sampling.

A second important detail is that the integer coordinates of the sample returned by
`PiecewiseConstant2D::Sample()` are used to index into `f` for
filter function evaluation. If
instead the point `p` was scaled up by the
size of the `f` array in each dimension and converted to an integer,
the result would occasionally differ from the integer coordinates computed during
sampling by `PiecewiseConstant2D` due to floating-point round-off error.
(Using the notation of Section 6.8.1, the issue is that with
floating-point arithmetic, .)
Again, variance would result, as the ratio might not be .

### 8.8.3 Box Filter

One of the most commonly used filters in graphics is the *box filter*
(and, in fact, when filtering and reconstruction are not addressed
explicitly, the box filter is the *de facto* result). The box filter
equally weights all samples within a square region of the image. Although
computationally efficient, it is just about the worst filter possible.
Recall from the discussion in Section 8.1.2 that the box
filter allows high-frequency sample data to leak into the reconstructed
values. This causes postaliasing—even if the original sample values were
at a high enough frequency to avoid aliasing, errors are introduced by poor
filtering.

Figure 8.51(a) shows a graph of the box filter, and Figure 8.52 shows the result of using the box filter to reconstruct two 1D functions.

For the step function we used previously to illustrate the Gibbs phenomenon, the box does reasonably well. However, the results are much worse for a sinusoidal function that has increasing frequency along the axis. Not only does the box filter do a poor job of reconstructing the function when the frequency is low, giving a discontinuous result even though the original function was smooth, but it also does an extremely poor job of reconstruction as the function’s frequency approaches and passes the Nyquist limit.

For this filter and all the following ones, we will not include the
rudimentary constructors and `Radius()` method implementations.

Evaluating the box filter requires checking that the given point is inside the box.

Sampling is also easy: the random sample `u` is used to linearly
interpolate within the filter’s extent. Since sampling is exact and the
filter function is positive, the weight is always 1.

Finally, the integral is equal to the filter’s area.

### 8.8.4 Triangle Filter

The triangle filter gives slightly better results than the box: the weight falls off linearly from the filter center over the square extent of the filter. See Figure 8.51(b) for a graph of the triangle filter.

Evaluating the triangle filter is simple: it is the product of two linear functions that go to 0 after the width of the filter in both the and directions. Here we have defined the filter to have a slope of , though the filter could alternatively have been defined to have a value of 1 at the origin and a slope that depends on the radius.

Because the filter is separable, its PDF is as well, and so each dimension
can be sampled independently. The sampling method uses a separate
`SampleTent()` utility function that is defined in
Section A.4.1.
Once again, the weight returned in
the `FilterSample` is always 1 because the filter is positive and
sampling is exact.

Finally, the triangle filter is easily integrated.

### 8.8.5 Gaussian Filter

Unlike the box and triangle filters, the Gaussian filter gives a reasonably
good result in practice. This filter applies a Gaussian bump that is
centered at the pixel and radially symmetric around it.
Figure 8.53 compares plots of the Gaussian
filter and the Mitchell filter (described in Section 8.8.6).
The Gaussian does tend to cause slight
blurring of the final image compared to some of the other filters, but this
blurring can help mask any remaining aliasing. This filter is the
default one used in `pbrt`.

The *Gaussian function* is parameterized by the position of the peak
and the standard deviation :

Larger values of cause a slower falloff, which leads to a blurrier image when the Gaussian is used as a filter.

The `GaussianFilter` is centered at the origin, so . Further,
the filter function subtracts the
value of the Gaussian at the end of its extent from the filter value in order
to make the filter go to 0 at its limit:

For efficiency, the constructor precomputes the constant term for in each direction.

The product of the two 1D Gaussian functions gives the overall filter value
according to Equation (8.26). The calls to
`std::max()` ensure that the value of 0 is returned for points outside
of the filter’s extent.

The integral of the Gaussian is

where is the error function. `GaussianIntegral()` evaluates its
value over a given range. The filter function’s integral can be
computed by evaluating the Gaussian’s integral over the filter’s
range and subtracting the integral of the offset that takes the filter to
zero at the end of its extent.

It is possible to sample from the Gaussian function using a polynomial
approximation to the inverse error function, though that is not sufficient
in this case, given the presence of the second term of the filter function in
Equation (8.26). `pbrt`’s
`GaussianFilter` implementation therefore uses a `FilterSampler`
for sampling.

### 8.8.6 Mitchell Filter

Filter design is notoriously difficult, mixing mathematical analysis and
perceptual experiments. Mitchell and Netravali (1988)
developed a family of parameterized filter functions in order to be able to
explore this space in a systematic manner. After analyzing test subjects’
subjective responses to images filtered with a variety of parameter values,
they developed a filter that tends to do a good job of trading off between
*ringing* (phantom edges next to actual edges in the image) and
*blurring* (excessively blurred results)—two common artifacts from
poor reconstruction filters.

Note from the graph in Figure 8.53(b) that
this filter function takes on negative values out by its edges; it has
*negative lobes*. In practice these negative regions improve the
sharpness of edges, giving crisper images (reduced blurring). If they
become too large, however, ringing tends to start to enter the image.
Furthermore, because the final pixel values can become negative, they
will eventually need to be clamped to a legal output range.

Figure 8.54 shows this filter reconstructing the two test functions. It does extremely well with both of them: there is minimal ringing with the step function, and it does a good job with the sinusoidal function, up until the point where the sampling rate is not sufficient to capture the function’s detail.

The Mitchell filter has two parameters called and . Although any values can be used for these parameters, Mitchell and Netravali recommend that they lie along the line .

The Mitchell–Netravali filter is the product of 1D filter functions in the and directions and is therefore separable.

The 1D function used in the Mitchell filter is an even function defined over the range . This function is made by joining a cubic polynomial defined over with another cubic polynomial defined over . This combined polynomial is also reflected around the plane to give the complete function. These polynomials are controlled by the and parameters and are chosen carefully to guarantee and continuity at , , and . The polynomials are

`Mitchell1D()` evaluates this
function. Its implementation is straightforward and is not included here.

As a cubic polynomial, sampling this filter function directly would require
inverting a quartic. Therefore, the `MitchellFilter` uses the
`FilterSampler` for sampling.

However, the function is easily integrated. The result is independent of the values of and .

### 8.8.7 Windowed Sinc Filter

Finally, the `LanczosSincFilter` class implements a filter
based on the sinc function. In practice, the sinc filter is often
multiplied by another function that goes to 0 after some distance. This
gives a filter function with finite extent. An additional parameter
controls how many cycles the sinc function passes through before it
is clamped to a value of 0.
Figure 8.55 shows a graph of three cycles of
the sinc function, along with a graph of the windowing function we
use, which was developed by Lanczos. The Lanczos window is
just the central lobe of the sinc function, scaled to cover the
cycles:

Figure 8.55 also shows the filter that we will
implement here, which is the product of the sinc function and the windowing
function. It is evaluated by the `WindowedSinc()` utility function.

Its implementation uses the `Sinc()` function, which in turn is
implemented using the numerically robust `SinXOverX()` function.

`LanczosSincFilter`.

Figure 8.56 shows the windowed sinc’s reconstruction results for uniform 1D samples. Thanks to the windowing, the reconstructed step function exhibits far less ringing than the reconstruction using the infinite-extent sinc function (compare to Figure 8.11). The windowed sinc filter also does extremely well at reconstructing the sinusoidal function until prealiasing begins.

The evaluation method is easily implemented in terms of the
`WindowedSinc()` function.

There is no convenient closed-form approach for sampling from the windowed
sinc function’s distribution, so a `FilterSampler` is used here as
well.

There is no closed-form expression of the filter’s integral, so its
`Integral()` method, not included in the text, approximates it using a
Riemann sum.