digital filtering. The largest body of antialiasing research stems from computer graphics
where high-quality rendering of complicated imagery is the central goal. The developed
algorithms have primarily addressed the tradeoff issues of accuracy versus efficiency.
Consequently, methods such as supersampling, adaptive sampling, stochastic sampling,
pyramids, and preintegrated tables have been introduced. These techniques are described
in Chapter 6.
-256 -192 -128
Figure 4.15: Bandlimited scanline appropriate for four-fold subsampling.
4.6 AIrlALIASING 111
-256 -192 -128 -64 0 64 128 192 256
l I I Till I [ [ [1 3 ß I Ill [I II
112 SAMPLING THEORY
This chapter has reviewed the basic principles of sampling theory. We have shown
that a continuous signal may be reconstructed from its samples if the signal is bandlim-
ited and the sampling frequency exceeds the Nyquist rate. These are the two necessary
conditions for image reconstruction to be possible. Since sampling can be shown to
replicate a signal's spectram across the frequency domain, ideal low-pass filtering was
introduced as a means of retaining the original spectrum while discarding its copies.
Unfortunately, the ideal low-pass filter in the spatial domain is an infinitely wide sine
function. Since this is difficult to work with, nonideal reconstruction filters are intro-
duced to approximate the reconstructed output. These filters are nonideal in the sense
that they do not completely attenuate the spectxa copies. Furthermore, they contribute to
some blurring of the original spectrum. In general, poor reconstruction leads to artifacts
such as jagged edges.
Aliasing refers to the phenomenon that occurs when a signal is undersampled. This
happens if the reconstruction conditions mentioned above are violated. In order to
resolve this problem, one of two actions may be taken. Either the signal can be bandlim-
ited to a range that complies with the sampling frequency, or the sampling frequency can
be increased. In practice, some combination of both options are taken, leaving some
relatively unobjectionable aliasing in the output.
Examples of the concepts discussed in this chapter are concisely depicted in Figs.
4.18 through 4.20. They attempt to illustrate the effects of sampling and low-pass filter-
ing on the quality of the reconstructed signal and its spectrum. The first row of Fig. 4.18
shows a signal and its spectxa, bandlimited to .5 cycle/pixel. For pedagogical purposes,
we txeat this signal as if it is continuous. In actuality, though, it is really a 256-sample
horizontal cross-section taken from the Mandrill image. Since each pixel has 4 samples
contributing to it, there is a maximum of two cycles per pixel. The horizontal axes of the
spectxa account for this fact.
The second row shows the effect of sampling the signal. Since fs = 1 sample/pixel,
there are four copies of the baseband speclrum in the range shown. Each copy is scaled
byfs= 1, leaving the magnitudes intact. In the third row, the 64 samples are shown con-
volred with a sine function in the spatial domain. This corresponds to a rectangular
pulse in the frequency domain. Since the sinc function is used here for image reconslruc-
tion, it must have an amplitude of unity value in order to interpolate the data. This forces
the height of the rectangular pulse in the frequency domain to vary in response to fs.
A few comments on the reciprocal relationship between the spatial and frequency
domains are in order here, particularly as they apply to the ideal low-pass filter. We
again refer to the variables A and W as defined in Section 4.3.3. As a sine function is
made broader, the value l/2W is made to change since W is decreasing to accommodate
zero crossings at larger intervals. Accordingly, broader sinc functions cause more blur-
ring and their spectxa reflect this by reducing the cut-off frequency to some smaller W.
Conversely, narrower sine functions cause less blurring and W takes on some larger
value. In either case, the amplitude of the sine function or its spectrum will change.
0 16 32 48 64
o 16 32 48 64 -2 -1 0 1 2
0 16 32 48 64 -2 -1 o 1 2
0 16 32 48 64 -2 -1 0 I 2
0 16 32 48 64 -2 -1 0 I 2
0 16 32 48 64 -2 -1 o I 2
Figure 4.18: Sampling and reconstruction (with an adequate sampling rate).
(Created by S. Feiner and G. Wolberg for [Foley 90]. Used with permission.)
That is, we can fix the amplitude of the sine function so that only the rectangular pulse of
the spectrum changes height A/2W as W varies. Altematively, we can fix A/2W to
remain constant as W changes, forcing us to vary A. The choice depends on the applica-
When the sine function is used to interpolate data, it is necessary to fix A to 1.
causing W to vary. This makes the amplitude of the spectrum's rectangular pulse change.
On the other hand, if the sine function is applied to bandlimit, not interpolate, the input
signal, then it is .important to fix A/2W to 1 so that the passband frequencies remain
intact. Since W is once again varying, A must change proportionately to keep A/2W con-
stant. Therefore, this application of the ideal low-pass filter requires the amplitude of the
sine function to be responsive to W.
In the examples presented below, our objective is to interpolate (reconstmc0 the
input and so A = 1 regardless of the sampling density. Consequently, the height of the
spectrum of the reconstruction filter changes. To make the Fourier transforms of the
filters easier to see, we have not drawn the frequency response of the reconstruction
filters to scale. Therefore, the rectangular pulse function in the third row of Fig. 14.18
actually has height A/2W= 1. The fourth row of the figure shows the result after apply-
ing the ideal low-pass filter. As sampling theory predicts, the output is identical to the
original signal. The last two rows of the figure illustrate the consequences of nonideal
reconstruction filtering. Instead of using a sine function, a triangle function correspond-
ing to linear interpolation was applied. In the frequency domain this corresponds to the
square of the sine function. Not surprisingly, the spectrum of the reconstructed signal
suffers in both the passband and the stopband.
The identical sequence of filtering operations is performed in Fig. 4.19. In this
figure, though, the sampling rate has been lowered to fs = .5, meaning that only one sam-
ple is collected for every two output pixels, Consequently, the replicated spectra are
multiplied by :5, leaving the magnitudes at 4. Unfortunately, this sampling rate causes
the replicated spectra to overlap. This, in turn, gives rise to aliasing, as depicted in the
fourth row of the figure. Applying the triangle function to perform linear interpolation
also yields poor results.
In order to combat these artifacts, the input signal must be bandlimited to accommo-
date the low sampling rate. This is shown in the second row of Fig. 14.20 where we see
that all frequencies beyond W=.25 are trancated. This causes the input signal to be
blurred. In this manner we have traded aliasing for blurring, a far less objectionable
artifact. Sampling this function no longer causes the replicated copies to overlap. Con-
volring with an ideal low-pass filter now properly isolates the bandlimited spectram.
0 16 32 48 64 -2 -1 0 1 2
0 16 32 48 64 -2 -1 0 1 2
0 16 32 48 64 -2 -i 0 I 2
0 16 32 48 64 .2 -I 0 I 2
Figure 4.19: Sampling and reconstruction (with an inadequate sampling rate).
(Created by S. Feiner and G. Wolberg for [Foley 90]. Used with permission.)
0 16 32 48 64
Image resampling is the process of txansforming a sampled image from one coordi-
nate system to another. The two coordinate systems are related to each other by the map-
ping function of the spatial transformation. This permits the output image to be gen-
erated by the following stxaightforward procedure. First, the inverse mapping function is
applied to the output sampling grid, projecting it onto the input. The result is a resam-
pling grid, specifying the locations at which the input is to be resampled. Then, the input
image is sampled at these points and the values are assigned to their respective output
The resampling process outlined above is hindered by one problem. The resam-
pling grid does not generally coincide with the input sampling grid, taken to be the
integer lattice. This is due to the fact that the range of the continuous mapping function
is the set of real numbers, a superset of the integer grid upon which the input is defined.
The solution therefore requires a match between the domain of the input and the range of
the mapping function. This can be achieved by convening the discrete image samples
into a continuous surface, a process known as image reconstruction. Once the input is
reconstxucted, it can be resampled at any position.
Conceptually, image resampling is comprised of two stages: image reconstruction
followed by sampling. Although resampling takes its name from the sampling stage,
image reconstruction is the implicit component in this procedure. It is achieved through
an interpolation procedure, and, in fact, the terms reconstruction and interpolation are
often used interehangeably.
The image resampling process is depicted in Fig. 5.1 for the 1-D case. A discrete
input (squares) is shown passing through the image reconstruction module, yielding a
continuous input signal (solid curve). Reconstruction is performed by convolving the
discrete input signal with a continuous interpolating function. The reconstructed input is
then modulated (multiplied) with a resampling grid (dashed arrows). Note that the
resampling grid is the result of projecting the output grid onto the input through a spatial
118 IMAGE REVAMPLING
Resamplingl i I -- I i I i
Grid ' '
I Spatial Transformation
Output ' ' i ' i
Grid I I i I I
ß ß ß ß
Input Samples Output Samples
transformation. After the reconstructed signal is sampled by the msampling grid, the
samples (cimles) are assigned to the uniformly spaced output image.
Image magnification and minification are two typical instances of image resam-
pling. These operations are known by many different names. For instance, stretching,
zooming, scaling up, interpolation, and upsampling are all informal terms used to
describe magnification. Similarly, minification, t compression, shrinking, scale reduc-
tion, decimation, and downsampling are all terms that describe the process of reducing
the size of an image. These two processes are illustrated in Fig. 5.2. In the top half of
the figure, the interval between two adjacent black and white pixels must be recon-
structed in order to generate five output points. A ramp is fitted between these points and
uniformly sampled at five locations to yield the intensity gradation appearing at the out-
put. In the bottom half of the figure, a scale reduction is shown. This was achieved by
discarding points, a method prone to aliasing. Later we shall review antialiasing algo-
rithms that use prefilters to bandlimit the input before resampling the continuous warped
signal. Prefilters will be shown to be related to the interpolation functions used in recon-
This term originated in the computer graphics literature [Smith 83].
Figure 5.2: Image magnification and minification.
form accurate image resampling. This chapter focuses on interpolation functions useful
in reconstructing a continuous function from sampled image data. Before proceeding to
image reconstruction, we briefly present an overview of ideal resampling. Although
somewhat theoretical, the presentation should serve to identify the roles of reconstruction
and prefiltering in their proper context. Together, they are used to define the ideal resam-
5.2. IDEAL IMAGE RESAMPLING
There are four basic elements to ideal image resampling: reconstraction, warping,
prefiltering, and sampling [Smith 83, Heckbert 89]. They are depicted in Fig. 5.3, and
outlined in Table 5.1.
The progression begins with f (u), the discrete input defined over integer values of
u. It is reconstructed into fc(U) through convolution with reconstruction filter r(u).
From sampling theory, we know that the ideal reconstruction filter is the sinc function.
The continuous input fc(U) is then warped according to mapping function m. The for-
ward map is given as x = rn (u) and the inverse map is u = m-1 (x). In this case, the warp
is defined as an inverse mapping. It is also possible to formulate this as a forward map-
ping instead. The spatial ramsformation leaves us with gc(x), the continuous warped
120 IMAGE RESAMPLING
f (a) g (x)
Discrete Input Discrete u x
l Reconstruct lSample
fc(U) gc(X) g(x)
Reconstructed InpUt u Warped Input TM x Continuous OutpUt x
Figure 5.3: Ideal resampling [Heckbert 89].
f(u), u Z
fc(U) = f(u)*r(u) = f(k)r(u-k)
gc(X) = fc(m-l(x))
g;(x) = go(x)* h(x) = f gc(t)h(x-t)dt
g(x) = g;(x)s(x)
Table 5,1: Elements of ideal resampling.
output. Depending on the inverse mapping function m-t(x), gc(X) may have arbitrarily
high frequencies. Therefore, it is bandlimited by function h (x) to conform to the Nyquist
rate of the output. The bandlimited result is g(x). This function is sampled by s (x), the
output sampling grid, to produce the discrete output g (x). Note that s (x), often referred
to as the comb function, is not required to sample the output at the same density as that of
There are only two filtering components to the entire resampling process: recon-
straction and prefiltering. We may cascade them into a single filter, derived as follows:
5.2 IDEAL IMAGE RE.SAMPLING 121
= ffc(m-(t)) h (x -t) dt
=llkzf(k)r(m-l(t)-k)] h(x-t) dt
= f(k) p(x,k)
p(x,k) = l r(m-l(t)-k) h(x-t) dt (5.2.2)
is the resarnpling filter that specifies the weight of the input sample at location k for an
output sample at location x [Heckbert 89].
Assuming that m (x) is invertible, we can express p(x,k) in terms of an integral in
the input space, rather than the output space. Substituting t = m (u), we have
p(x,k) = lr(u-k)h(x-m(u)) u du (5.2.3)
where I m/Ou I is the determinant of the Jacobian matrix interrelating the input and out-
put coordinate systems. In one dimension,
In two dimensions,
Ca= ;; (5.2.5)
where xu = Ox/Ou, and similar notation holds for the other partial derivatives.
Either the input-space or output-space integral can be used to define the resampling
filter. In the input-space form, p is expressed in terms of a reconstruction filter and a
warped prefilter. This can be readily justified by noting that the reconstruction filter is
applied before the warp and therefore it can be applied directly to the input. The
prefilter, however, is applied after the warp and so its domain, still defined in terms of u,
must undergo the geometric transformation. Since equal increments in u do not generally
correspond to identical increments in re(u), the prefilter is warped. This formulation of
the resampling filter is depicted in Fig. 5.4. A similar ease holds for the output-space
form of p, which is written in terms of a warped reconstruction filter and a prefilter.
Therefore, the actual warping is incorporated into either the reconstruction filter or
prefilter, but not both.
The resampling filter takes on a simple form for space-invariant linear warps. In
that case, the resampling filter can be shown to be equivalent to the convolution of the
reconstruction filter and prefilter [Heckbert 89]. Expressed in input-space form, we have
122 IMAGE RE,qAMPLING
f(R ....... tion] fc(u).-----l] g(u) g(x)
Figure 5.4: Ideal resampling with input-space resampling filter [Heckbert 89].
p(x,k) = p'(m-(x)-k) (5.2.6)
= h'(u)* r(u)
= Jlh(uJ)] *r(u)
where J is the Jacobian matrix and u = m-l(x ) -k. This formulation is suitable for linear
warps defined in terms of forward mapping functions, i.e., m (u) = uJ.
In the special case of magnification, we may ignore the profilter altogether, treating
it instead as an impulse function. This is due to the fact that no high frequencies are
introduced into the output upon magnification. Conversely, minification introduces high
frequencies and does not require any reconstruction of the input image. Consequently,
we can ignore the reconstruction filter and treat it simply as an impulse function. There-
Pmax(x,k) = r(m l(x)-k) (5.2.7a)
p,,in(X,/C) = I JIh (x -m )) (5.2.7b)
Equations (5.2.7a) and (5.2.7b) lead us to an important observation about the shape
of reconstruction filters and profilters for linear warps. According to Eq. (5.2.7a), the
shape of the reconstraction filter does not change in response to the mapping function.
Indeed, magnification is achieved by selecting a reconstruction filter and directly con-
volring it across the input. Its shape remains fixed independently Of the magnification
scale factor. A similar procedure is taken in minification, whereby a reconstruction filter
is replaced by a prefilter. The prefilter is selected on the basis of some desired frequency
response characteristics. Unlike reconstruction filters, though, the actual shape mast be
scaled by an amount linearly related to the minification factor. As the input is increas-
ingly decimated, the prefilter must become broader and shorter. It becomes broader in
order to average more neighboring pixels together, thereby further bandlimiting the
input. Since larger neighborhoods are used to compute each output pixel, the normalized
5.2 IDEAL IMAGE RESAMPLING 123
weights applied to the input decrease to reflect the diminishing impact of each input sam-
pie. As a result, the prefilter grows shorter.
This observation is a direct consequence of the reciprocal relationship between the
spatial and frequency domains. Due to the importance of this property, a proof is
presented below. We start by writing the expression for the Fourier transform of h (u).
h(u) f h(u)e-i2nfUdu (5.2.8)
Note that we use the symbol here to denote a transform pair. After we warp the
input h (u) through mapping function m (u), we get
h (m (u)) I h (m (u)) e 12ndu (5.2.9)
Letting x = au = m (u) and dx = u du, we have
h(m(u)) f h(x)e -i2npn (x) dx (5.2.10)
where m - (x) = x/a and ]m/u I = IJ I = a. This gives us
h (au) > 1 f h (x) e-i2nfx/adx (5.2.11)
This equation expresses the reciprocal relationship between the spatial and frequency
domains. Notice that multiplying the spatial axis by a factor of a results in dividing the
frequency axis and the spectrum values by that same factor.
This proves to be a fundamental result in linear filtering theory that bears significant
consequences. For instance, we would ideally like to use narrow filters in the spatial
domain. In this manner, each output pixel can be computed by weighting only a small
number of input samples. However, the reciprocal relationship tells us that narrow filters
in the spatial domain correspond to wide frequency spectrums. This, however, is
undesirable as it hinders our attempts to avoid aliasing due to spectral overlaps. On the
other hand, wide spatial filters are costly, but they do permit as to perform more effective
bandlimiting. This tradeoff between narrow filters in the spatial domain and good filter
response in the frequency domain is at the heart of filter design.
The remainder of this chapter focuses on interpolation for reconstruction, a central
component of image resampling. This area has received extensive treatment due to its
practical significance in numerous applications. Although theoretical limits on image
reconstruction are derived by sampling theory, the algorithms proposed in this chapter
address tradeoff issues in accuracy and complexity.
124 IMAGE RESAMPLING
Interpolation is the process of determining the values of a function at positions lying
between its samples. It achieves this process by fitting a continuous function through the