Digital image warping

f (x)

fo f l f

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

4 ,_

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

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

k=j-2

yields

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

f2 141

(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 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) (

(C) (d)

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-06

(C) (d)

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)

.7s ...h....;..-..L....h., ,..i.......i......L....i-.

IHCf)}

0.0011

le05

-4-3-2-101234 -4-3-2-101234

½)

Figure 5.13: (a) Hamming window; (b) Windowed sinc; (c) Spectrum; (d) Log plot.

The Blackman window is similar to the Hann and Hamming windows. It is defined

as

i 2rx 4rx N - 1

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.01

1-04 [

le-06 L

(C) (d)

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

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

le-04]

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

0.1

0.e01

lo-05

-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

h4(x)

0.01

lc-04 J

tL

%)

o.oot I

re-06

(C) 

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