IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. I8, NO. II , NOVEMBER 1999
`
`!049
`
`Survey: Interpolation Methods
`1n Medical Image Processing
`
`Thomas M. Lehmann,* Member, IEEE, Claudia Gonner, and Klaus Spitzer
`
`Abstract- Image interpolation techniques often are required
`in medical imaging for image generation (e.g., discrete back
`projection for inverse Radon transform) and processing such as
`compression or resampling. Since the ideal interpolation function
`spatially is unlimited, several interpolation kernels of finite size
`have been introduced. This paper compares 1) truncated and win(cid:173)
`dowed sine; 2) nearest neighbor; 3) linear; 4) quadratic; S) cubic
`B-spline; 6) cubic; g) Lagrange; and 7) Gaussian interpolation
`and approximation techniques with kernel sizes from 1 x 1 up to
`8 x 8. The comparison is done by: 1) spatial and Fourier analyses;
`2) computational complexity as well as runtime evaluations; and
`3) qualitative and quantitative interpolation error determinations
`for particular interpolation tasks which were taken from common
`situations in medical image processing.
`For local and Fourier analyses, a standardized notation is intro(cid:173)
`duced and fundamental properties of interpolators are derived.
`Successful methods should be direct current (DC)-constant and
`interpolators rather than DC-inconstant or approximators. Each
`method's parameters are tuned with respect to those properties.
`This results in three novel kernels, which are introduced in this
`paper and proven to be within the best choices for medical
`image interpolation: the 6 x 6 Blackman-Harris windowed sine
`interpolator, and the CZ-continuous cubic kernels with N = 6
`and N = 8 supporting points.
`For quantitative error evaluations, a set of SO direct digital
`X rays was used. They have been selected arbitrarily from
`clinical routine. In general, large kernel sizes were found to be
`superior to small interpolation masks. Except for truncated sine
`interpolators, all kernels with N = 6 or larger sizes perform
`significantly better than N = 2 or N = 3 point methods
`(p ~ 0.005). However, the differences within the group of large(cid:173)
`sized kernels were not significant. Summarizing the results, the
`cubic 6 x 6 interpolator with continuous second derivatives,
`as defined in (24), can be recommended for most common
`interpolation tasks. It appears to be the fastest six-point kernel
`to implement computationally. It provides eminent local and
`Fourier properties, is easy to implement, and has only small
`errors. The same characteristics apply to B-spline interpolation,
`but the 6 x 6 cubic avoids the intrinsic border effects produced
`by the B-spline technique.
`However, the goal of this study was not to determine an
`overall best method, but to present a comprehensive catalogue of
`methods in a uniform terminology, to define general properties
`and requirements of local techniques, and to enable the reader
`to select that method which is optimal for his specific application
`in medical imaging.
`
`Manuscript received August 6, 1998; revised September 24, 1999. This
`work was supported in part by the German Research Community (DFG) under
`Grants Re 427/5-1, Sp 538/1-2, and Schm 1268/1-2. The Associate Editor
`responsible for coordinating the review of this paper and recommending its
`publication was A. Manduca. Asterisk indicates corresponding author.
`*T. M. Lehmann, C. Ganner, and K. Spitzer are with the Institute of Medical
`Informatics, Aachen University of Technology (RWTH), 52057 Aachen,
`Germany (e-mail: lehmann@computer.org).
`Publisher Item Identifier S 0278-0062(99)10280-5.
`
`Index Terms- Approximation, B-splines, cubic polynomials,
`image resampling, interpolation.
`
`I. INTRODUCTION
`
`I MAGE interpolation has many applications in computer
`
`vision. It is the first of the two basic resampling steps
`and transforms a discrete matrix into a continuous image.
`Subsequent sampling of this intermediate result produces
`the resampled discrete image. Resampling is required for
`discrete image manipulations, such as geometric alignment and
`registration, to improve image quality on display devices or
`in the field of lossy image compression wherein some pixels
`or some frames are discarded during the encoding process
`and must be regenerated from the remaining information
`for decoding. Therefore, image interpolation methods have
`occupied a peculiar position in medical image processing [1].
`They are required for image generation as well as in image
`post-processing. In computed tomography (CT) or magnetic
`resonance imaging (MRI), image reconstruction requires in(cid:173)
`terpolation to approximate the discrete functions to be back
`projected for inverse Radon transform. In modem X-ray imag(cid:173)
`ing systems such as digital subtraction angiography (DSA),
`interpolation is used to enable the computer-assisted alignment
`of the current radiograph and the mask image. Moreover,
`zooming or rotating medical images after their acquisition
`often is used in diagnosis and treatment, and interpolation
`methods are incorporated into systems for computer aided di(cid:173)
`agnosis (CAD), computer assisted surgery (CAS), and picture
`archieving and communication systems (PACS).
`Image interpolation methods are as old as computer graphics
`and image processing. In the early years, simple algorithms,
`such as nearest neighbor or linear interpolation, were used for
`resampling. As a result of information theory introduced by
`Shannon in the late 1940's, the sine function was accepted
`as the interpolation function of choice. However, this ideal
`interpolator has an infinite impulse response (IIR) and is not
`suitable for local interpolation with finite impulse response
`(FIR). From the mathematical point of view, Taylor or La(cid:173)
`grange polynomials have been suggested to approximate the
`sine function [2]. This is documented in most textbooks on
`numerical analysis [3]. Thereafter, due to their numerical
`efficiency, different families of spline functions have been
`used instead.
`A great variety of methods with confusing naming can be
`found in the literature of the 1970' s and 1980's. B-splines
`sometimes are referred to as cubic splines [4], while cubic
`
`0278--0062/99$10.00 © 1999 IEEE
`
`0001
`
`exocad GmbH, et al.
`Exhibit 1010
`
`Authorized licensed use limited to: TEL AVIV UNIVERSITY. Downloaded on January 21, 2010 al 02:38 from IEEE Xplore. Restrictions apply.
`
`

`

`1050
`
`IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 18, NO. 11, NOVEMBER 1999
`
`TABLE I
`PREVIOUS PAPERS COMPARING MORE THAN THREE INTERPOLATION METHODS
`
`Interpolation sheme
`
`Truncated sine
`Windowed sine
`Nearest neighbor
`Linear
`Quadratic (approx.)
`Quadratic (interpo.)
`B-spline (approx.)
`13-spline (interpol.)
`Cubic, 2 x 2
`Cubic, 4 x 4
`Cubic, 6 x 6
`Cubic, 8 x 8
`Lagrange
`Gaussian
`
`[4]
`Ac
`
`Ac
`Ac
`
`[5]
`
`[6]
`
`[7]
`
`[8]
`
`[10]
`
`AB
`ABc
`
`Ac
`AB
`AB A Be ace A Be
`
`A Be A
`abc
`
`AB
`
`AB
`ace
`A
`ABcd AB ABc ace
`ad
`
`A Be
`ABcd
`
`[11]
`
`[12]
`
`(13]
`AB
`AC
`ABC A Be
`ABc ABC ABc
`A Be ABC A Be
`ABc
`ABc
`ABc
`
`A
`ac
`
`[14]
`
`acD
`aBcCD
`
`ad
`aBcCD
`
`A
`ABc AC
`
`ABc
`
`aBcCD
`
`ace
`
`ABC
`
`Abbreviations (a) kernels' derivation, (A) including plots; (b) Fourier analyP1i11, (B) including plots; (c) image based qualitative comparison
`by 8ubjects and (C) quantitative interpolation crmr determination; (d) complexity eva.luation and (D) runtime measurements. Note that this
`paper covers the cubic 2-point and 8-poinl methods with 'abcCdD, and all other rnws with 'ABcCdD' .
`
`the computation time of cubic kernels to 60% by the use of
`interpolation is also known as cubic convolution [5], [6],
`quadratic functions yielding similar quality [11].
`high-resolution spline interpolation [7], and bi-cubic spline
`Table I summarizes previous work comparing interpolation
`interpolation [8], [9]. In 1983, Parker, Kenyon, and Troxel
`methods. In addition to Appledorn's and Dodgson's recent pro(cid:173)
`published the first paper entitled "Comparison of Interpo(cid:173)
`posals, most comparative studies include neither the windowed
`lation Methods" [7], followed by a similar study presented
`sine technique nor the Lagrange method and also exclude
`by Mealand in 1988 [6]. However, previous work of Hou
`large kernels for cubic interpolation with 6 x 6 or 8 x
`and Andrews, as well as that of Keys, also compare global
`8 supporting points. We will see below that those methods
`and local interpolation methods ([ 4] and [5], respectively)
`(Table I). The Fourier transform was used in these studies to
`perform superiorly in most applications.
`In medical diagnostic applications, not only the kernel's
`evaluate different 2 x 2 and 4 x 4 interpolation methods.
`frequency properties must be taken into account but also
`Parker et al. pointed out that, at the expense of some increase
`the appearance of images after resampling. Keep in mind
`in computing time, the quality of resampled images can be
`that many imaging systems violate the sampling theorem and
`improved using cubic interpolation when compared to nearest
`introduce aliasing. Unser, Aldroubi, and Eden asked subjects
`neighbor, linear, or B-spline interpolation. However, to avoid
`to rank as lifelike the magnified Lena test image in descending
`further perpetuation of misconceptions, which have appeared
`repeatedly in the literature, it might be better to refer to
`order [10]. Although this type of evaluation particularly de(cid:173)
`pends on the images and their geometric transforms, visually
`their B-spline technique as B-spline approximation instead
`assessed interpolation quality was found to be important for
`of interpolation. Maeland named the correct (natural) spline
`interpolation as B-spline interpolation and found this technique
`kernel selection [4], [5], [7], [11], [13], [14], [16], [20]. Others
`use the Fourier power spectrum of their test images before and
`to be superior to cubic interpolation [6].
`after interpolation to determine the quality of their technique
`In more recent reports, not only hardware implementations
`[8], [17]. Alternatively, Schaum suggests an error spectrum for
`for linear interpolation [15] and fast algorithms for B-spline
`interpolation [10] or special geometric transforms [8], [9], [14],
`comparing the performance of various interpolation methods
`[12]. Table I also summarizes characteristics used by previous
`but also nonlinear and adaptive algorithms for image zoom(cid:173)
`authors to compare interpolation methods.
`ing with perceptual edge enhancement [16], [17] have been
`However, other bases for evaluating the appropriateness of a
`published. However, smoothing effects are most bothersome
`given interpolation scheme may be indicated, depending on the
`if large magnifications are required. In addition, shape-based
`task. For example, the visual performance can be quantified
`and object-based methods have been established in medicine
`in terms of similarity or sharpness [17], [21]. Furthermore,
`for slice interpolation of three-dimensional (3-D) data sets
`magnifications by factors of eight [10] or even two seldom
`[18]. In 1996, Appledorn presented a new approach to the
`are required in many clinical applications and, thus, comparing
`interpolation of sampled data [ 19]. His interpolation functions
`this capability may be irrelevant if not inappropriate for cer(cid:173)
`are generated from a linear sum of a Gaussian function and
`tain tasks. Accordingly, determining the optimal interpolation
`their even derivatives. Kernel sizes of 8 x 8 were suggested.
`method for a variety of specific clinical applications is still a
`Again, Fourier analysis was used to optimize the parameters
`problem.
`of the interpolation kernels. Contrary to large kernel sizes
`This paper presents a comprehensive survey of existing
`and complex interpolation families causing a high amount
`image-interpolation methods. They are expressed using a stan(cid:173)
`of computation, the use of quadratic polynomials on small
`dardized terminology and are compared by means of lo-
`regions was recommended by Dodgson in 1997. He reduces
`0002
`
`Authorized licensed use limited to: TEL AVIV UNIVERSITY. Downloaded on January 21, 2010 al 02:38 from IEEE Xplore. Restrictions apply.
`
`

`

`LEHMANN et al.: SURVEY: INTERPOLATION METHODS
`
`1051
`
`cal and Fourier analyses, qualitative and quantitative error
`determinations, computational complexity evaluations, and
`run time measurements. Three representative interpolation
`tasks have been selected from clinical routine for comparing
`performance. These are described in detail in the next section.
`Relative performance and task-specific dependencies likewise
`are examined with regard to variations between image types
`and transform parameters. The interpolation methods as well
`as their parameters and variations are presented and discussed
`in Section III. The results of analyzing Fourier properties,
`interpolation errors, and run time are presented in Section IV
`and discussed in the last section of this paper.
`
`II. INTERPOLATION TASKS IN MEDICAL IMAGING
`
`Image resampling is required for every geometric transform
`of discrete images except shifts over integer distances or
`rotations about multiples of 90°. Geometric transforms differ
`with respect to their complexity. Usually, affine transforms
`such as magnifications [4], [7], [10], [11], [16], [17], [20], [22]
`or rotations [l], [8], [9], [14] are used to evaluate interpolation
`methods. High enlargements up to 4; 5.25; and 8 times ([11],
`(16], [22]; [5]; and [4], (10], [22], respectively) as well as 10
`or 16 repetitions of the transform ((20] and [14] respectively)
`are used to enhance the blurring effects incorporated with in(cid:173)
`terpolation. In this paper, simple expansions in one dimension,
`different rotations, and complex perspective transforms within
`the range of clinical applications are used to compare the inter(cid:173)
`polation methods. Representative clinical images produced by
`a variety of diagnostic modalities provide a reasonable basis
`for evaluation. These include medical photographs, magnetic
`resonance displays, and digital radiographs.
`
`A. Correction of Aspect Ratios in CCD-Photographs
`Magnifications in only one dimension are required to correct
`aspect ratios of digital photographs acquired with CCD-sensors
`and frame grabber cards. To preserve the original information,
`one dimension of the image is expanded, rather than shrinking
`the other.
`Fig. 1 shows a digital photograph of a human eye. The
`positions of the Purkinje reflections within the pupil are used
`for strabometry [23]. The picture was selected because of its
`uniform histogram and its excellent sharpness. The defraction
`patterns of the reflections in the pupil and the eyelashes are
`particularly sharp and thus provide an interpolation challenge.
`The correction of its 4/3-aspect ratio requires the expansion
`of the x axis by 1.3. Note that in this task the pixels in every
`third column of the expanded image can be taken from the
`original data without any modification.
`
`(a)
`
`(b)
`
`Fig. I. CCD-arrays or frame grabber often require the correction of their
`aspect ratio. (a) In this example, the image of a human eye was acquired for
`strabometl)' [23]. (b) The 4/3 expansion in x direction was performed by
`linear interpolation.
`
`after resampling. In performing rotations on a discrete grid,
`neither image rows nor columns can be reproduced unless the
`rotation angle is a multiple of 90°. Therefore, this task results
`in larger interpolation errors than those produced by a simple
`correction of aspect ratio.
`
`C. Perspective Projection of X-Ray Images
`In intraoral radiology, the geometric transform of one ra(cid:173)
`diograph into another from the same dental region acquired at
`a different time is described by a perspective projection (25].
`Each pixel (x, y) within the reference image is transformed
`into the position ( x', y') in the subsequent image
`a.1x + a2y + a3
`/ ~x + ar,y + au
`d
`,
`B. Rotation of MRI Sections
`x=
`an y=
`a7x + asy + 1
`a1x + asy + 1
`Fig. 2 shows an MRI of the human head. In multimodal
`image registration, CT or MRI slices must be transformed
`Suppose Fig. 3(a) was acquired with an imaging plate
`as an affine projection to fit data from functional imaging
`perpendicular to the beam cone and the viewpoint is behind the
`modalities, such as positron emission tomography (PET) or
`image plate looking toward the X-ray tube. Fig. 3(b) illustrates
`single photon emission computed tomography (SPECT) [24].
`the rotation of the image plate with its upper right comer
`By restricting affine transformations to rotations, the number of
`moving toward the X-ray tube. Although distortions in clinical
`image points remains approximately constant both before and
`routine are less drastic, the parameters ai of the perspective
`0003
`
`(1)
`
`Authorized licensed use limited lo: TEL AVIV UNIVERSITY. Downloaded on January 21, 2010 at 02:38 from IEEE Xplore. Restrictions apply.
`
`

`

`1052
`
`IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 18, NO. 11, NOVEMBER 1999
`
`(a)
`
`(b)
`
`Interpolation is required to rotate discrete images. The MR image (a)
`Fig. 2.
`was 45° rotated using linear interpolation (b).
`
`projection have been chosen such that the number of pixels
`containing information is roughly halved after the projection.
`Therefore, the interpolation errors will be greatest within the
`examples presented here.
`
`lll. INTERPOLATION METHODS
`
`For image resampling, the interpolation step must recon(cid:173)
`struct a two-dimensional (2-D) continuous signal s(x, y) from
`its discrete samples s(k, l) with s, x, y E lR and k, l E
`JN°. Thus, the amplitude at the position (x, y) must be
`estimated from its discrete neighbors. This can be described
`formally as the convolution of the discrete image samples
`with the continuous 2-D impulse response 2oh(x, y) of a 2-D
`reconstruction filter
`
`s(x, y) = EE s(k, l) · wh(x - k, y - l).
`
`(2)
`
`k
`l
`Usually, symmetrical and separable interpolation kernels are
`used to reduce the computational complexity
`wh(x, y) = h(x) · h(y).
`
`(a)
`
`(b)
`
`Intraoral X rays are acquired with respect to perspective projection.
`Fig. 3.
`Therefo1 e, perspective tra.nsfonns must be pcrfonncd for automatical adjust(cid:173)
`ment. The parameters a; in (1) are chosen as follows: a 1 = 1.25, a 2 = 0.35,
`a3 = 1.10, a1 = 0.20, as = 0.80, as = 5.00, a1 = 0.002, as = -0.0006.
`This equals a movement of the upper right comer of the image plate toward
`the X-ray tube. The 340 x 256 original pixels (a) are reduced to only 39 770
`pixels (45.7%) which still contain image information after interpolation (b).
`
`Fig. 4. One-dimensional decomposition of the 2-D N x N interpolation of
`the point (x, y) . (Reprinted with permission from [l].)
`
`first. The small grey intermediate points in Fig. 4 are generated
`by four one-dimensional (1-0) interpolations. They are used
`for the final 1-D interpolation in the y direction.
`
`(3) A. Ideal Interpolation
`Fig. 4 illustrates the interpolation of the point (x,y) in a 4 x
`Following the sampling theory, the scanning of a continuous
`image s(x, y) yields infinite repetitions of its continuous
`4 neighborhood. Interpolation is performed in the x direction
`0004
`
`Authorized licensed use limited to: TEL AVIV UNIVERSITY. Downloaded on January 21, 2010 at 02:38 from IEEE Xplore. Restrictions apply.
`
`

`

`LEHMANN el al: SURVEY: INTERPOLATION METHODS
`
`1053
`
`• 4
`
`. 2
`
`• 8
`
`• 6
`
`.4
`
`.2
`
`1.
`
`0.1
`
`0.01
`
`0.001
`
`0 . 0001
`
`(b)
`(a)
`Ideal interpolation. (a) Kernel plotted for /.-.: / < 3. (b) Magnitude of Fourier transfonn. (c) Logarithmic plot of magnitude.
`
`(c)
`
`Fig. 5.
`
`- 0
`
`-~
`
`fO
`
`- 10
`
`- 5
`
`5
`
`10
`
`.8
`
`.6
`
`. 1
`
`.2
`
`(b)
`
`1.
`
`O. l
`
`0.01
`
`0.001
`
`0.0001
`
`-1 0
`
`-5
`
`6
`
`10
`
`(c)
`
`(a)
`
`Fig. 6. Truncated sine interpolation, N = 5. (a) Kernel. (b) Magnitude of Fourier transform. (c) Logarithmic plot of magnitude.
`
`spectrum S( u, v) in the Fourier domain, which do not overlap
`since the Nyquist criterion is satisfied. If this is so, and
`only then, the original image s(x, y) can be reconstructed
`perfectly from its samples s(k, I) by multiplication of an
`appropriate rectangular prism in the Fourier domain. The 1-
`D ideal interpolation equals the multiplication with a rect
`function in the Fourier domain and can be realized in the
`spatial domain by a convolution with the sine function
`
`Some fundamental properties of any interpolator can be
`derived from this ideal interpolation function. Idcalh(x) is
`positive from zero to one, negative from one to two, positive
`
`guarantee that the image is not modified if it is resampled on
`the same grid. Therefore, kernels satisfying
`
`from two to three, and so on. For h(O) = 1 these zero crossings
`h(O) = 1
`
`{
`
`h(x) = 0,
`
`Jxl = 1, 2, . ..
`
`(5)
`
`I<lealh( ) -
`( )
`.
`sin( 7rX) -
`x - - --
`-smcx.
`1rX
`
`(4)
`
`Fig. 5(a) shows the ideal IIR-interpolator I<lealh(x ). The plot
`was truncated within the interval -3 < x < 3. The magnitude
`JI<leal H(w )J of the Fourier transform I<leal H(w) of the infinite
`kernel Ideal h( x) is plotted within the interval -47r :S w =
`27rf::::; 47r is shown in Fig. 5(b). The interval -11" < w < 7r
`is called passband and f = 1/2 or w = 7r the cutoff
`point or Nyquist frequency. The transfer function of the ideal
`interpolator is constant and one in the passband. In addition, a
`logarithmical plot of the filter's Fourier response is presented
`in Fig. 5(c) to emphasize ripples in the stopband Jwl > 7r.
`The ideal kernel's transfer function is zero valued within
`the stopband. Note that Figs. 5-24 are constructed in similar
`fashion using the same scaling and plotting conventions.
`
`avoid smoothing and preserve high frequencies. They are
`called interpolators. We will see below that better suited
`kernel functions tend to have this general shape. In contrast
`to interpolators, kernels that do not fulfill (5) are named
`approximators. Note that this strict distinction is not always
`reflected in literature.
`Sampling the interpolated (continuous) image is equivalent
`to interpolating the (discrete) image with a sampled interpo(cid:173)
`lation function [7]. The sampling of the interpolation function
`aliases the higher frequencies of the interpolation function into
`the lower ones. In the case of ideal interpolation only, higher
`frequencies do not exist and therefore, within the interval
`-0.5 < f < 0.5, the sampled interpolation function has the
`same Fourier spectrum as the unsampled function. However, it
`is necessary to examine not only the continuous interpolation
`0005
`
`Authorized licensed use limited lo: TEL AVIV UNIVERSITY. Downloaded on January 21, 2010 at 02:38 from IEEE Xplore. Restrictions apply.
`
`

`

`1054
`
`IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 18, NO. II, NOVEMBER 1999
`
`3
`
`. 6
`
`-~
`
`.2
`
`l(l
`
`1.
`
`0 . 1
`
`0.01
`
`0.001
`
`0.0001
`
`~-~
`i! :
`: . L ... ___.__._._._._.
`
`-10
`
`-5
`
`0
`
`10
`
`:1
`
`(b)
`(a)
`Fig. 7. Truncated sine interpolation, N = 6. (a) Kernel. (b) Magnitude of Fourier transform. (c) Logarithmic plot of magnitude.
`
`(c)
`
`0 6
`
`0 . 4
`
`0 8
`
`. 6
`
`.4
`
`0.2
`
`l.
`
`0 . 1
`
`0 . 01
`
`0.001
`
`0.0001
`
`- 3
`
`-
`
`2
`
`3
`
`-0.2
`
`- 0
`
`5
`
`5 ··To
`
`(a)
`(b)
`(c)
`Fig. 8. Blackman-Harris windowed sine interpolation N = 6. (a) Kernel. (b) Magnitude of Fourier transform. (c) Logarithmic plot of magnitude.
`
`function h(x) but also typically sampled interpolation func- where
`tions h(k). Particularly, the sum of all samples should be one
`for any displacement 0 :=::; d < 1
`00 L h(d+k)=1.
`
`(6)
`
`k=-=
`
`III(x) = L 8(x + k)
`
`00
`
`k=-=
`
`and, by definition, the weight of a single delta impulse
`corresponds to the amplitude of the kernel h at the position
`of 8. Recognizing the Russian letter "scha," the train of delta
`functions is named the scha-function III(x). If the integrand in
`1r f:n, which equals 1 for .f = 0, one can
`(7) is extended by e-i2
`discover the definition of the Fourier transform in (7). Then,
`the function to be Fourier transformed is h( x) ·III( x + d) and
`from ( 6) and (7) we obtain
`H(f) ><III(f)·ei21rfdl = 1 ¢= {H(O) =. 1
`H(J) = o,
`' -
`J=O -
`
`Ill= 1,2, ...
`(8)
`where H(J) denotes the Fourier transform of h(x). Because
`the conditions in (8) are not sufficient but necessary in the
`context of interpolation, they are used to distinguish DC(cid:173)
`constant from DC-inconstant kernels in the Fourier domain.
`
`This means that for any displacement d the direct current (DC)(cid:173)
`amplification will be unity and the energy of the resampled
`image remains unchanged. In other words, the mean brightness
`of the image is not affected if the image is interpolated
`or resampled. Therefore, kernel functions that satisfy or fail
`condition (6) are named DC-constant or DC-inconstant, re(cid:173)
`spectively. The next sections will show that superior kernels
`are DC-constant.
`Equation (6) also is called the partition of unity condition
`(20], which easily can be evaluated in the Fourier domain.
`Referring to information theory, the sum (6) of discrete sam(cid:173)
`ples of the kernel h(x) equals the area under the continuous
`function obtained by multiplying (or sampling) h(x) with a
`train of delta functions 8(x)
`
`f h(d + k) = 100
`
`k=-oo
`
`h(x) · III(x + d)dx
`-oo
`
`(7)
`
`B. Sine Interpolation
`Although the sine function provides an exact reconstruction
`of s(x, y), it spatially is unlimited. There are two common
`0006
`
`Authorized licensed use limited to: TEL AVIV UNIVERSITY. Downloaded on January 21, 2010 at 02:36 from IEEE Xplore. Restrictions apply.
`
`

`

`LEHMANN et al.: SURVEY: INTERPOLATION METHODS
`
`1055
`
`0 8
`
`1
`
`2
`
`3
`
`-0.2
`
`(a)
`
`1.
`
`0.1
`
`J.01
`
`. 001
`
`)001
`
`(b)
`
`(c)
`
`-lO -5
`
`C
`
`~
`
`10
`
`Fig. 9. Nearest neighbor interpolation. (a) Kernel. (b) Magnitude of Fourier transform. (c) Logarithmic plot of magnitude.
`
`1
`
`1.
`
`0.1
`
`0.01
`
`0.001
`
`0.0001
`
`-3
`
`- 2
`
`2
`
`3
`
`-0.2
`
`(a)
`
`(b)
`
`0
`
`s
`
`0
`
`(c)
`
`Fig. JO. Linear interpolation. (a) Kernel. (b) Magnitude of Fourier transform. (c) Logarithmic plot of magnitude.
`
`approaches for overcoming this drawback, truncation and
`windowing with a window function w(x) = const(x) = 1
`and w(x) -:f. const(x), respectively
`
`0:::; lxl < N/2
`elsewhere
`
`(9)
`
`within the passband. With respect to the passband properties
`of a truncated kernel, odd numbers of supporting points are
`preferable.
`Another idea to make the sine function usable for spatial
`convolution might be to use it with a less severe window
`w(x) than the rect function. Ostuni et al. discuss the use
`of a cosine function w(x) = cos('rrx/N) for reslicing fMRI
`data [26]; Schaum uses a Hanning window, which is just a
`raised cosine, to taper the interpolation kernel's edges and
`remove Gibbs's overshoot in the transform [12]; and Walberg
`compares several window functions for interpolation with
`windowed sine kernels by Fourier analysis [13]. A systematic
`approach on the use of windows for harmonic analysis with
`the discrete Fourier transform is given by Harris who declared
`the Kaiser-Bessel and Blackman-Harris windows to be the top
`performers [27]. When using the three-term Blackman-Harris
`window
`
`( 4x)
`211"" N + w2 co.s
`
`211"" N
`
`(10)
`
`( 2x)
`
`( ) - { Idealh(x). w(x),
`Sinch
`O
`-
`N X
`'
`where N denotes the number of the finite kernel's supporting
`points. By definition, SinchN(x) fulfills the requirement (5). In
`other words, all windowed or truncated sine kernels necessarily
`are real interpolators.
`Truncation is equivalent to the multiplication of Idealh(x)
`with a rectangular function in the spatial domain, which is
`tantamount to a convolution with a sine function in the fre(cid:173)
`quency domain. Therefore, truncations of the ideal interpolator
`produce ringing effects in the frequency domain because a
`considerable amount of energy is discarded. Figs. 6 and 7
`demonstrate this effect, which also is referred to as the Gibbs's
`phenomenon [4], produced by a truncated sine function with
`N = 5 and N = 6 supporting points, respectively. In addition,
`the partition of unity condition (8) is violated by any choice
`of N < oo. In other words, all truncated sine kernels are DC(cid:173)
`inconstant. The area of the function differs more from one for
`even kernel sizes than for odd. Therefore, raising the kernel
`from N = 5 to N = 6 significantly enlarges the overshoots
`0007
`
`w(x) = w0 + w1 cos
`
`with N = 6 and
`
`Wo = 0.42323
`W1=0.49755
`W2 = 0.07922
`
`Authorized licensed use limited to: TEL AVIV UNIVERSITY. Downloaded on January 21, 2010 at 02:38 from IEEE Xplore. Restrictions apply.
`
`

`

`1056
`
`IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 18, NO. 11, NOVEMBER 1999
`
`1
`
`0 . 8
`
`0 6
`
`.~
`
`- 3
`
`- 2
`
`3
`
`-0.2
`
`- 0
`
`-5
`
`10
`
`1.
`
`0. 1
`
`0. 01
`
`0.001
`
`0 . 0001
`
`·~------J
`
`0
`
`10
`
`- 10
`
`- 5
`
`(b)
`(a)
`Fig. 11 . Quadratic approximation, a = 1/2. (a) Kernel. (b) Magnitude of Fourier transform. (c) Logarithmic plot of magnitude.
`
`(c)
`
`0 6
`
`0 4
`
`-3
`
`- ?.
`
`-0.2
`
`1.
`
`O. l
`
`0.001
`
`0.0 0 01
`
`- 10
`
`- 5
`
`0
`
`s 1 0
`
`(a)
`Fig. 12. Quadratic interpolation, a = 1. (a) Kernel. (b) Magnitude of Fourier transform. (c) Logarithmic plot of magnitude.
`
`(b)
`
`(c)
`
`a DC-constant interpolator is obtained. Note that most other
`window functions, including those used in [12] and [13] result
`in kernels that do not have this superior property.
`Fig. 8 shows the Blackman-Harris windowed sine kernel.
`The kernel ' s half wave between 2 S lxl < 3 is suppressed
`significantly in comparison with the ideal or truncated kernels.
`Therefore, the ripples in the stopband are below 0.01 %, but
`higher frequencies within the passband are attenuated also.
`The largest gain within the stop band is 0.5 at the cutoff point.
`
`C. Nearest Neighbor Interpolation
`The easiest way to approximate the sine function by a
`spatially limited kernel is given by the nearest neighbor
`method. The value s( x) at the location ( x) is chosen as the
`next known value s(k). Therefore, only N = 1 supporting
`point is required for the nearest neighbor interpolation. This
`is tantamount to convolution with a rect function [Fig. 9(a)]
`0 S lxl < 0.5
`elsewhere.
`
`(11)
`
`sidelobes in those regions of the frequency domain where the
`repetitions of S caused by scanning s should be suppressed
`[Fig. 9(c)]. The gain in the passband rapidly falls off to
`2/'rr :::::! 64% at the cutoff point, and the amplitude of the
`side maxima is more than 20%. Therefore, strong aliasing
`and blurring effects are associated with the nearest neighbor
`method for image interpolation.
`
`D. Linear Interpolation
`For separated bi-linear interpolation, the values of both
`direct neighbors are weighted by their distance to the opposite
`point of interpolation. Therefore, the linear approximation of
`the sine function follows the triangular function
`
`h (x) = {1 - Ix! , 0 S lxl < 1
`
`0,
`
`elsewhere.
`
`2
`
`(12)
`
`The triangular function h2 ( x) corresponds to a modest
`low-pass filter H 2 (f) in the frequency domain (Fig. 10).
`Again, h2(0) = 1, h2(±1, 2, ... ) = 0, H2 (0) = 1, and
`H 2(±1, 2, . .. ) = 0. Therefore, the linear kernel is a DC(cid:173)
`Clearly, h1 (x) is a DC-constant interpolator.
`constant interpolator. The sidelobes in the stopband are below
`Fig. 9(b) shows that the Fourier spectrum of the nearest
`10%, which still is considerable. Therefore, the main disad(cid:173)
`neighbor kernel equals the sine function (expressed in the
`frequency domain). The logarithmical scale shows prominent
`vantages of linear interpolation are both the attenuation of
`0008
`
`Authorized licensed use limited lo: TEL AVIV UNIVERSITY. Downloaded on January 21, 2010 at 02:38 from IEEE Xplore, Restrictions apply.
`
`

`

`LEHMANN et al: SURVEY: INTERPOLATION METHODS
`
`1057
`
`1;
`
`0.8
`
`- 3
`
`-2
`
`2
`
`3
`
`-0.2
`
`(a)
`
`0.
`
`0 6
`
`0 4
`
`5
`
`10
`
`(b)
`
`0.1
`
`0.01
`
`0.001
`
`0.0001
`
`- 10
`
`, __
`
`_,
`
`0
`
`(c)
`
`\ !\
`''l
`
`ll_
`
`5
`
`.10
`
`Fig. 13. Cubic B-spline approximation. (a) Kernel. (b) Magnitude of Fourier transform. (c) Logarithmic plot of magnitude.
`
`high-frequency components and the aliasing of the data beyond
`the cutoff point into the low frequencies [7].
`
`a DC-constant kernel. Hence, the following four equations are
`required to establish appropriate values for the five remaining
`parameters:
`
`E. Quadratic Approximation
`One of the most frequently applied concepts to create sine
`like interpolation kernels is the use of algebraic polynomials.
`Their advantage is easy determination and uniform approxima(cid:173)
`tion of continuous functions at finite intervals. In the previous
`sections, constant and linear pol