f (x)
fo f l f
Xo X 1 X 2 X 3 X 4 X 5 X 6
Figure 5.9: A spline consisting of 6 piecewise cubic polynomials.
The four coefficients offt can be defined in terms of the data points and their first
(or second) derivatives. Assuming that the data samples are on the integer lattice, each
134 IMAGE RESAMPLING
spaced one unit apart, then the coefficients, defined in terms of the data samples and their
first derivatives, are given below.
a0 = Y (5.4.16a)
a = y (5.4.16b)
a2 = 3Ay, - 2y -Y,+I (5.4.16c)
a3 = -2Ay +y +Y+I (5.4.16d)
where Ay = y+ - y.
Although the derivatives are not supplied with the data, they are derived by solving
the following system of linear equations.
41
141
y
y
1 4 _
4 ,_
-Sy0 + 4yl +Y2
30'2 -Y0)
30'3 -Yl)
30'n-t - Yn-3)
-Y-3 - 4yn-2 + 5yn-
(5.4.17)
The not-a-knot boundary condition [de Boor 78] was used above, as reflected in the
first and last rows of the matrices. It is superior to the artificial boundary conditions com-
monly reported in the literature, such as the natural or cyclic end conditions, which have
no relevance in our application. Note that the need to solve a linear system of equations
arises from global dependencies introduced by the constxaints for continuous first and
second derivatives at the knots. A complete derivation is given in Appendix 2.
In order to compare interpolating cubic splines with other methods, we must
analyze the interpolation kernel. Thus far, however, the piecewise interpolating polyno-
mials have been derived without any reference to an interpolation kernel. We seek to
express the interpolating cubic spline as a convolution in a manner similar to the previous
algorithms. This can be done with the use of cubic B-splines as interpolation kernels
[Hou 78].
5.4.5.1. B-Splines
A B-spline of degree n is derived through n convolutions of the box filter, B 0.
Thus, B t =B0*B0 denotes a B-spline of degree 1, yielding the familiar triangle filter
shown in Fig. 5.7a. Interpolation by B consists of a sequence of stxaight lines joined at
the knots continuously. This is equivalent to linear interpolation.
S.4 INTERPOLATION KERNELS 13
The second-degree B-spline B2 is produced by convolving Bo*B 1. Using B2 to
interpolate data yields a sequence of parabolas that join at the knots continuously
together with their slopes. The span orB2 is limited to three points.
The cubic B-spline B 3 is generated from convolving Bo*B2. That is,
B 3 = Bo*Bo*Bo*B o. The interpolation with B 3 is composed of a series of cubic poly-
nomials that join at the knots continuously together with their slopes and curvatures, i.e.,
their first and second derivatives. Figure 5.10 summarizes the shapes of these low-order
B-splines.
-l.5 -.5 .5 1.5
-2 -1 0 1 2
Figure 5.10: Low-order B-splines are derived from repeated box filters.
Denoting the cubic B-spline interpolation kernel as h, we have the following piece-
wise cubic polynomials defining the kemel.
[31xl3-6lxl 2+4 0-< Ixl < l
h(x) = -}lx13+61x12-121xl+8 l
2_< Ixl
This kemel is sometimes called the Parzen window.
There are several properties of cubic B-splines worth noting. As in the cubic con-
volution method, the extent of the cubic B-spline is over four points. This allows two
points on each side of the centxal interpolated region to be used in the convolution. Con-
sequently, the cubic B-spline is shift-invariant as well.
Unlike cubic convolution, however, the cubic B-spline kernel is not interpolatory
since it does not satisfy the necessary consmint that h (0)= 1 and h(1)= h(2)= 0.
Instead, it is an approximating function that passes near the points but not necessarily
through them. This is due to the fact that the kernel is strictly positive.
136 IMAGE RESAMPLING
The posifivity of the cubic B-spline kernel is actually attractive for our image pro-
cessing application. When using kernels with negative lobes, (e.g., the cubic convolution
and windowed sinc functions), it is possible to generate negative values while interpolat-
ing positive data. Since negative intensity values are meaningless for display, it is desir-
able to use strictly positive interpolation kernels to guarantee the positivity of the interpo-
lated image.
There are problems, however, in dkectly interpolating the data with kernel h, as
given in Eq. (5.4.18). Due to the low-pass (blur) characteristics of h, the image under-
goes considerable smoothing. This is evident by examining its frequency response where
the stopband is effectively suppressed at the expense of additional attenuation in the
passband. This leads us to the development of an interpolation method built upon the
local support of the cubic B-spline.
5.4.5.2. Interpolating B-Splines
Interpolating with cubic B-splines requires that at data point x/, we again satisfy Eq.
(5.4.7). Namely,
j+2
f(xj) = ch(xj-x,) (5.4.19)
k=j-2
From Eq. (5.4.18), we have h(0)=4/6, h(-1)= h(1)= 1/6, and h(-2)= h(2)= 0. This
yields
1
f (xj) = (cj_ 1 q- 4cj + Cj+l) (5.4.20)
Since this must be true for all data points, we have a chain of global dependencies for the
ck coefficients. The resulting linear system of equations is similar to that obtained for the
derivatives of the cubic interpolating spline algorithm. We thus have,
f0 [4 1
fl i41
f2 141
C2
(5.4.21)
Labeling the three matrices above as F, K, and C, respectively, we have
F = K C (5.4.22)
The coefficients in C may be evaluated by multiplying the known data points F with the
inve[se of the tridiagonal matrix K.
C = K - F (5.4.23)
5.4 INTERPOLATION KERNELS 137
The inversion of tridiagonal matrix K has an efficient'algorithm that is solvable in
linear time [Press 88]. In [Lee 83], the matrix inversion step is modified to introduce
high-frequency emphasis. This serves to compensate for the undesirable low-pass filter
imposed by the point-spread function of the imaging system.
In all the previous methods, the coefficients c were taken to be the data samples
themselves. In the cubic spline interpolation algorithm, however, the coefficients must
be determined by solving a tridiagonal matrix problem. After the interpolation
coefficients have been computed, cubic spline interpolation has the same computational
cost as cubic convolution.
5.4.6. Windowed Sinc Function
Sampling theory establishes that the sine function is the ideal interpolation kernel.
Although this interpolation filter is exact, it is not practical since it is an IIR filter defined
by a slowly converging infinite sum. Nevertheless, it is perfectly reasonable to consider
the effects of using a trancated, and therefore finite, sinc function as the interpolation ker-
nel.
The results of this operation are predicted by sampling theory, which demonstrates
that huncation in one domain leads to ringing in the other domain. This is due to the fact
that truncating a signal is equivalent to multiplying it with a rectangle function Rect(x),
defined as
Rect(x) = .5 -< Ix l (5.4.24)
Since multiplication in one domain is convolution in the other, lynncation amounts to
convolving the signal's spectram with a sinc function, the transform pair ofRect (x). We
have already seen an example of this in Fig. 4.7. Since the stopband is no longer elim-
inated, but rather attenuated by a ringing filter (i.e., a sinc), the input is not bandlimited
and aliasing artifacts are introduced. The most typical problems occur at step edges,
where the Gibbs phenomena becomes noticeable in the form of undershoots, overshoots,
and ringing in the vicinity of edges. In [Ratzel 80], the author found this method to per-
form poorly.
The Rect function above served as a window, or kemel, that weighs the input signal.
In Fig. 5.11a, we see the Rect window extended over three pixels on each side of its
center, i.e., Rect(6x) is plotted. The corresponding windowed sinc function h(x) is
shown in Fig. 5.1lb. This is simply the product of the sine function with the window
function, i.e., sinc(x)Rect(6x). Its spectrum, shown in Fig. 5.11c, is nearly an ideal
low-pass filter. Although it has a fairly sharp mmsifion from the passband to the stop-
band, it is plagued by ringing. In order to more clearly see the values in the spectrum, we
use a logarithmic scale for the vertical axis of the spectram in Fig. 5.11 d. The next few
figures will be illustrated by using this same four-part format.
Ringing can be mitigated by using a different windowing function exhibiting
smoother fall-off than the rectangle. The resulting windowed sine function can yield
138 IMAGE REVAMPLING S.4 INTERPOLATION KERNELS 139
Rect (x)
-.2s L--.i..-...!-....--i...-.-i.-.-..i.-.-..!...-...::......i.....-:: ]
4-3-2-101234 4-33-101234
(a) (
-4 -3 -2 -1 0 1 2 3 4 -4 -3 -2 -1 0 1 2 3 4
(C) (d)
Figure 5.11: (a) Rectangular window; (b) Windowed sinc; (c) Spectrum; (d) Log plot.
better results. However, since slow fall-off requires larger windows, the computation
remains costly.
Aside from the rectangular window mentioned above, the most frequently used win-
dow functions are: Harm, t Hamming, Blackman, and Kaiser [Antoniou 79]. These filters
identify a quantity known as the ripple ratio, defined as the ratio of the maximum side-
lobe amplitu.d.e to the main-lobe amplitude. Good filters will have small ripple ratios to
achieve effective attenuation in the stopband. A tradeoff exists, however, between ripple
ratio and main-lobe width. Therefore, as the ripple ratio is decreased, the main-lobe
width is increased. This is consistent with the reciprocal relationship between the spatial
and frequency domains, i.e., narrow bandwidths correspond to wide spatial functions.
In general, though, each of these smooth window functions is defined over a small
finite extent. This is tantamount to multiplying the smooth window with a rectangle
function, While this is better than the Rect function alone, there will inevitably be some
form of aliasing. Nevertheless, the window functions described below offer a good
compromise between tinging and blurring.
' Due to Julius yon Harm. It is often mistakenly referred to as the Htinning window.
5.4.6.1. Hann and Hamming Windows
The Hann and Hamming windows are defined as
{ 2for N- 1
+(1-a)cos2 i- Ixl < 2
HannlHamming(x) = - (5.4.25)
otherwise
where N is the number of samples in the windowing function. The two windowing func-
tions differ in the choice of x. In the Hann window x=0.5, and in the Hamming window
0=0.54. Since they both amount to a scaled and shifted cosine function, they are also
known as the raised cosine window.
The spectra for the Hann and Hamming windows can be shown to be the sum of a
sinc, the spectrum of Rect(x), with two shifted counterparts: a sinc shifted to the tight by
2x/(N - 1), as well as one shifted to the left by the same amount. This serves to cancel
the right and left side lobes in the specmm of Rect(x). As a result, the Hann and Ham-
ming windows have reduced side lobes in their spectra as compared to those of the rec-
tangular window. The Hann window is illustrated in Fig. 5.12.
-4 -3 -2 -I 0 1 2 3 4 -4 -3 -2 -1 0 1 2 3 4
(a) (b)
IH(f)l IHCf)l
0.0011
le-04 I
le-05
le-06
-4-3-2-1 0 1 2 3 4 -4-3-2-1 0 1 2 3 4
(C) (d)
Figure 5.12: (a) Hann window; (b) Windowed sinc; (c) Spectram; (d) Log plot.
IIq] I
140 IMAGE RESAMPL1NG
Notice that the passband is only slightly attenuated, but the stopband continues to retain
high frequency components in the stopband, albeit less than that of Rect(x). It performs
somewhat better in the stopband than the Hamming window, as shown in Fig. 5.13. This
is partially due to the fact that the Hamming window is discontinuous at its ends, giving
rise to "kinks" in the spectrum.
-4 -3 -2 -1 0 1 2 3 4 -4 -3 -2 -I 0 1 2 3 4
(a) (b)
IHCf)[
.7s ...h....;..-..L....h., ,..i.......i......L....i-.
IHCf)}
0.0011
le-04 I
le05
le-06
-4-3-2-101234 -4-3-2-101234
½)
Figure 5.13: (a) Hamming window; (b) Windowed sinc; (c) Spectrum; (d) Log plot.
5.4.6.2. Blackman Window
The Blackman window is similar to the Hann and Hamming windows. It is defined
as
i 2rx 4rx N - 1
ß 42+0.5cos----T+0.08cos---ZT Ix{ < 2
Blackman (x) = - - (5.4.26)
otherwise
The purpose of the additional cosine term is to further reduce the ripple ratio. This win-
dow function is shown in Fig. 5.14.
5.4 INTERPOLATION KERNELS 141
Blackman (x )
-4 -3 -2 -1 0 1 2 3 4 -4 -3 -2 -1 0 1 2 3 4
(a) (b)
IHtf)l
0.1
0.01
o.ol ]
1-04 [
1½-0 ]
le-06 L
-4-3-2-1 0 1 2 3 4 -4-3-2-1 0 1 2 3 4
(C) (d)
Figure 5.14: (a) Blackman window; (b) Windowed sinc; (c) Spectrum; (d) Log plot.
5.4.6.3. Kaiser Window
The Kaiser window is defined as
f 1o(b) N- 1
Kaiser(x) =l: (0 'xl< 2 (5.4.27)
otherwise
where ot is a free parameter and
f f 211t'
I0 is the zeroth-order Bessel function of the first kind. This can be evaluated to any
desired degree of accuracy by using the rapidly converging series
Io(n) = 1 + (5.4.29)
k=l
142 IMAGE RESAMPLING
The Kaiser window leaves the filter designer much flexibility in controlling the rip-
ple ratio by adjusting the parameter at. As at is incremented, the level of sophistication of
the window function grows as well. Therefore, the rectangular window corresponds to a
Kaiser window with at = 0, while more sophisticated windows such as the Hamming win-
dow correspond to o = 5. This formulation facilitates a tradeoff between ringing and
edge softening.
5.4.6.4. Lanczos Window
Windowed sinc functions are notorious for producing ringing artifacts near edges.
Although they are an improvement over truncated sinc functions, they retain a fairly
sharp transition from passband to stopband. Superior filters can be designed by imposing
further constraints on the filter response in the frequency domain.
Reasonable constraints to impose on the kernel include: unity gain in the low-pass
region with cut-off at frequency fi, zero gain at high frequencies beyond f2, and linear
fall-off in the transition range between fi and f2. This frequency response can be
expressed as the convolution of two boxes. In the spatial domain, this corresponds to the
multiplication of two sinc functions, yielding a function known as the Lanczos window.
The widths of the two sinc functions determine the extent of the transition range.
The two-lobed Lanczos window function is defined as
sin(;c/2) 0 -< Ix I < 2
Lanczos 2(x) = :x/2 (5.4.30)
0 2-
The Lanczos2 window function is the central lobe of a sinc function. It is wide
enough to extend over two lobes of the ideal low-pass filter, i.e., a second sinc function.
The windowed sinc function is therefore given by the product sinc (x)Latczos2(x). This
can be rewritten as sinc (x) sinc (x/2) Rect(x/4), where the first term is the ideal low-pass
filter, the second term is Lanczos 2(x), and Rect(x/4) is the rectangular function that tan-
cares Lanczos2 past x=2. Note that its abscissa is x/4 because Rect is defined over
-.5 < x < .5. The spectrum of this product is Rect(f)*Rect(2f)*sinc (4f), where * is
convolution. The Lanczos 2(x) window function is shown in Fig. 5.15.
This formulation can be generalized to an N-lobed window function by replacing
the value 2 in Eq. (5.4.30) to the value N. For instance, the 3-lobed Lanczos window is
defined as
sin(mr/3) 0 -< I x I < 3
Lanczos3(x) = rx/3 (5.4.31)
0 3-
The Lanczos3(x) window function is shown in Fig. 5.16. As we let more more lobes
pass under the Lanczos window, then the spectrum of the windowed sinc function
becomes Rect(f)*Rect(Nf)*sinc(2Nf). This proves to be a superior frequency
5.4 INTERPOLATION KERNELS 143
Irtf)l
0.1
0.01 [
le-04]
le,-05 [
-4 -3 -2 -1 0 1 2 3 4 -4 -3 -2 -1 0 1 2 3 4
(C) (d)
Figure 5.15: (a) Lanczos2 window; (b) Windowed sinc; (c) Spectrum; (d) Log plot.
response than that of the 2-lobed Lanczos window because the Rect(Nf) term causes
faster fall-off in the transition region and sinc (2N f) is a narrower sinc function that pro-
duces less deleterious ringing artifacts.
5.4.6.5. Gaussian Window
The Gaussian function is defined as
1 e_X2/2o2 (5.4.32)
Gauss(x) = 2q'o
where o is the standard deviation. Ganssians have the nice property that their spectrum
is also a Gaussian. They can be used to directly smooth the input for prefihering pur-
poses or to smooth the sinc function for windowed sinc reconstruction filters. Further-
more, since the tails of a Gaussian diminish rapidly, they may be truncated and still pro-
duce results that are not plagued by excessive ringing. The rate of fall-off is determined
by o, with low values of resulting in faster decay.
The general form of the Gaussian may be expressed more conveniently as
Gausso(x) = 2 -(x)2 (5.4.33)
144 IMAGE RESAMPLING
-4 -3 -2 -1 0 1 2 3 4 -4 -3 -2 -1 0 1 2 3 4
(a) (b)
IHf)l
0.1
0.01
0.e01
1½-04
lo-05
lo-06
-4-3-2-1 0 1 2 3 4 -4-3-2-1 0 1 2 3 4
(C) (d)
Figure 5.16: (a) Lanczos3 window; (b) Windowed sinc; (c) Spectram; (d) Log plot.
Two plots of this function are shown in Fig. 5.17 with = 1/2 and 1/'. The latter is a
wider Gaussian than the first, and its magnitude doesn't become negligible until two sam-
ples away from the center. The benefit of these choices of ct is that many of their
coefficients for commonly used resampling ratios are scaled powers of two, which makes
way for fast computation. In [Turkowskl 88a], these two functions have been examined
for use in image resampling.
-2 -1 0 1 2
Figure 5.17: Two Gaussian functions.
5.4 INTERPOLATION KERNELS 145
5.4.7. Exponential Filters
A superior class of reconstruction filters can be derived using exponential functions.
Consider, for instance, the hyberbolic tangent function tanh defined in Eq. (5.4.34).
e x -- e-x
tanh (x) = -- (5.4.34)
e x + e -x
This function has several desirable properties. First, it converges quickly to + 1. Second,
its transition from -1 to 1 is sharp. We can sharpen the transition even further by scaling
the domain, i.e., use tanh(kx) for k Z 1. In addition, this function is infinitely differenti-
able everywhere, i.e., it satisfies an important smoothness constraint. These properties
are readily apparent in Fig. 5.18, which illustrates tanh (kx) for k = 1, 4, and 10. Notice
that the function quickly approximates Rect for larger values of k.
-2 -1 0 1 2
Figure 5.18: Scaled hyperbolic tangent function.
Given tanh (kx) as our starting point, we can define a new function that resembles
the ideal low-pass filter Rect(f), i.e., a box in the frequency domain. This is done by
treating tanh as one half of Rect, and then merely compositing that with a mirror image
of itself. Since tanh lies between -1 and 1, some care must be taken to normalize the
expression so that it yields a box of unity height. The resulting function is given as
Hk(f)=[tanh(k(5+fc))+l'l[tanh(k(-+fc))+l- 1 (5.4.35)
where fc is the cut-off frequency. In our examples, we shall use f½ = .5 to conform to the
Nyquist rate. The purpose of the addition and division operations is to normalize H,(f)
so that 0 < Hk(f) < 1.
Function H,t is treated as the desired spectrum of our reconstruction filter. By vary-
ing k, we can control the shape of the spectram. For low values of k, H is smooth and
resembles a Gaussian function. As k is made larger, Hk will have increasingly sharper
comers, eventually approximating a Rect function. Figure 5.19 shows Hk(f ) for k = 1, 4,
and 10.
Having established H, to be the desired spectrum of our interpolation kemel, the
actual kemel is derived by computing the inverse Fourier transform of Eq. (5.4.35). This
146 IMAGE nESAMPLING
-2 -1 0 1 2
Figure 5.19: Spectrum HkO e ) is a function of tanh (kx).
gives us hk(x), as shown in Fig. 5.20. Not surprisingly, it has infinite extent. However,
unlike the sine function that decays at a rate of l/x, hk decays exponentially fast. This is
readily verified by inspecting the log plots in Fig. 5.20? This means that we may truncate
it with negligible penalty in reconstruction quality. The truncation is, in effect, implicit
in the decay of the filter. In practice, a 7-point kernel (with 3 points on each side of the
center) yields excellent results [Massalin 90].
n4(x)
-4 -3 -2 -1 0 1 2 3 4
(a)
h4(x)
0.1
0.01
0.001 [
lc-04 J
le-05 I
tL
-4 -3 -2 -1 0 1 2 3 4
%)
h10(x) h10(x)
o.oot I
le-05
re-06
-4-34-101234 -4-34-101234
(C)
Figure 5.20: Interpolation kernels derived from H4(f) and Ht0Oe).
t Note that a linear fall-off in log scale conesponds to an exponential function.
COMPARISON OF INTERPOLATION METHODS 147
5.5. COMPARISON OF INTERPOLATION METHODS
The quality of the popular interpolation kernels are ranked in ascending order as fol-
lows: nearest neighbor, linear, cubic convolution, cubic spline, and sine function. These
interpolation metheds are compared in many sources, including [Andrews 76, Parker 83,
Maeland 88, Ward 89]. Below we give some examples of these techniques for the
magnification of the Star and Madonna images. The Star image helps show the response
**Share with your friends:** |