This work presents a new method of fast cubic and higher order image interpolation. The evaluation of the piecewise n-th order polynomial kernels is accelerated by transforming the polynomials into the interval [0,1], which has the advantage that some terms of the polynomials disappear, and that several coefficients could be precalculated, which is proven in the paper. The results are exactly the same as using standard n-th order interpolation, but the computational complexity is reduced. Calculating the interpolation weights for the cubic convolution only needs about 60% of the time compared to the classical method optimized by the Horner's rule. This allows a new efficient implementation for image interpolation

Bischof, Horst

Perko, Roland

DSpace at University of West Bohemia

EFFICIENT IMPLEMENTATION OF HIGHER ORDERIMAGE INTERPOLATIONRoland PerkoGraz University of TechnologyInffeldgasse 16A–8010 GrazAustriaperko@icg.tu-graz.ac.atHorst BischofUniversity of Technology GrazInffeldgasse 16A–8010 GrazAustriabischof@icg.tu-graz.ac.atABSTRACTThis work presents a new method of fast cubic and higher order image interpolation. The evaluation of thepiecewise n-th order polynomial kernels is accelerated by transforming the polynomials into the interval[0, 1], which has the advantage that some terms of the polynomials disappear, and that several coefficientscould be precalculated, which is proven in the paper. The results are exactly the same as using standardn-th order interpolation, but the computational complexity is reduced. Calculating the interpolationweights for the cubic convolution only needs about 60% of the time compared to the classical methodoptimized by the Horner’s rule. This allows a new efficient implementation for image interpolation.KeywordsInterpolation, Reconstruction, Performance1 INTRODUCTIONImage interpolation is the process of reconstruct-ing a spatially continuous image from a set of dis-crete equidistant samples. It is fundamental inmany image processing operations, such as mag-nification, subpixel translation, rotation, deforma-tion or warping. These general operations requireimage values at locations from which no sample isavailable. From sampling theory it is well knownthat the sinc-function is the ideal interpolationkernel, which, however, cannot be used in prac-tise [Unser99]. In order to obtain acceptable re-construction in terms of computational speed andmathematical precision, it is required to design akernel that is of finite extent and approximatesthe sinc-function as much as possible.This is commonly done by symmetrical piecewisen-th order polynomial kernels, where the classi-cal cubic convolution kernel was first introducedPermission to make digital or hard copies of all or part ofthis work for personal or classroom use is granted withoutfee provided that copies are not made or distributed forprofit or commercial advantage and that copies bear thisnotice and the full citation on the first page. To copy oth-erwise, or republish, to post on servers or to redistributeto lists, requires prior specific permission and/or a fee.WSCG SHORT Communication papers proceedingsWSCG’2004, February 2–6, 2004, Plzen, Czech Republic.Copyright UNION Agency – Science Pressto image processing by [Rifma73]. A quantitativeevaluation of these interpolation kernels can befound in [Meije99], [Meije01] and [Lehma99].For arbitrary image warping every pixel of the re-sulting image has to be interpolated, which is ofhigher complexity than the evaluation of the geo-metric warping function itself. Therefore a speedup of the image interpolation task is desired andpresented in this paper.The paper is structured as follows. First, the con-cept of 1-D signal reconstruction using symmet-ric piecewise n-th order polynomial kernels is pre-sented (section 2). Next, the transformation of thepolynomials is described (section 3). In section4 an analysis of theoretical and practical perfor-mance of the proposed method is given. Finally,concluding remarks are made in section 5.2 RECONSTRUCTIONAs shown in [Meije99] the symmetrical piecewisen-th order polynomial kernel, which approximatesthe sinc function in [−(n− 1), (n− 1)], is definedas follows:h(x) ={∑nj=0 aij |x|j if i ≤ |x| < i+ 10 otherwise(1)To insure continuous, smooth interpolation, it isnecessary to set constraints at the knots. Flatfield interpolation requires that∑∞−∞ h(x−m) =1, where m = (n + 1)/2. Which leads to twoconstraints:h(0) = 1 and h(x) = 0 for |x| = 1, · · · ,m− 1 (2)h(l)(x) must be continuous at |x| = 0, 1, · · · ,m(3)With l = 0, 1, · · · , k where k = 0 for n = 1 andk = n− 2 for n > 1 which is proven in [Meije99].2.1 Cubic kernel exampleThe cubic 1-D kernel function is build up ofthird-order polynomials and approximates thesinc-function in the interval [−2, 2]. The kernel isgiven as:hC(x) =a03|x|3 + a02|x|2 + a01|x|+ a00 if 0 ≤ |x| < 1a13|x|3 + a12|x|2 + a11|x|+ a10 if 1 ≤ |x| < 20 otherwise (4)With the defined constraints in equation 2 and3 we get seven equations in eight unknowns. Byallowing a13 = α to be a tuneable parameter, thesystem can be solved, yielding following values:a03 = (α + 2), a02 = −(α + 3), a01 = 0, a00 =1, a13 = α, a12 = −5α, a11 = 8α and a10 =−4α. The free parameter α is chosen that way,that it satisfies the flatness constraint which leadsto α = − 12 [Keys81]. The cubic convolution kernelby Keys is given by the following equation and isshown in figure 1.p0(x) =32|x|3 − 52|x|2 + 1 if 0 ≤ |x| < 1p1(x) = −12 |x|3 +52|x|2 − 4|x|+ 2 if 1 ≤ |x| < 2(5)The standard way to calculate the interpolationweights is to evaluate the piecewise polynomialsfunctions on the supporting points. Every poly-nomial is evaluated for the according distance ξ.In case of cubic convolution these distances are ξ,1− ξ, ξ + 1 and 2− ξ where ξ ∈ [0, 1].3 TRANSFORMATION OFTHE POLYNOMIALSTo avoid the calculation of these distances all poly-nomials are transformed in the space [0, 1]. This isdone by evaluating the polynomial p0(x) at ξ = xand ξ = 1 − x, p1(x) at ξ = x + 1 and ξ = 2 − xanalytically for the cubic case. For higher (n-th)order kernels the strategy is the same and we getn+ 1 polynomials, for i odd and k = i/2:fi(x) = pk(x+ k)fi+1(x) = pk((k + 1)− x) (6)Withfi(x) =n∑j=0aijxj x ∈ [0, 1] (7)The cubic convolution kernel is shown in figure1 and figure 2 shows the transformed cubicpolynomials.−2 −1 0 1 200.10.20.30.40.50.60.70.80.91p1(2−x)p0(1−x)p0(x)p1(x+1)Figure 1: Key’s cubic convolution kernel.In order to save computational time we can proofthe following 5 lemmas about the properties of thetransformed polynomial coefficients. The proofsare given in the appendix. Figure 3 shows thestatements of the lemma graphically on the exam-ple of quintic interpolation coefficients.Lemma 1: The lowest order coefficients ai0 are0 for all polynomials fi(x), except for the firstpolynomial f0(x) where the coefficient is 1.Lemma 2: The highest order coefficients ainof pairwise symmetrical polynomials, fi(x) andfi+1(x), i odd, have the same values with oppositesign.0 0.5 100.10.20.30.40.50.60.70.80.91f0(x)0 0.5 100.10.20.30.40.50.60.70.80.91f1(x)0 0.5 100.10.20.30.40.50.60.70.80.91f2(x)0 0.5 100.10.20.30.40.50.60.70.80.91f3(x)Figure 2: Key’s cubic convolution kernel trans-formed to [0, 1].Lemma 3: All coefficients of order j ≤ n − 2of pairwise non symmetrical polynomials havethe same coefficients of odd order and the samecoefficients with opposite sign for even order.Lemma 4: All coefficients anj with j ≤ n − 2 ofthe polynomial fn(x) are 0.Lemma 5: All coefficients a0j with odd or-der j with j ≤ n − 2 are 0 for the polynomialf0(x).Figure 3: Statements of the given lemma shown onthe example of quintic (5-th order) interpolationcoefficients.Using these 5 lemmas we can formulate follow-ing efficient calculation. The distance x has thesame value for all polynomials, so the powers ofx can be precalculated. Furthermore the multipli-cations with a power of x and the pairwise samecoefficients have to be calculated only once andcoefficients with zeros could be skipped.3.1 Example for the cubic caseThe classical cubic convolution kernel accordingto Keys is given by equation 5.Transforming the polynomials to the space [0, 1]gives 4 polynomials:f0(x) =32x3 − 52x2 + 1f1(x) = −32x3 + 2x2 +12xf2(x) = −12x3 + x2 − 12xf3(x) =12x3 − 12x2 (8)According to the observations of the transformedpolynomials, here following values are calculated:x2, x3, 32x3, − 12x3 and 12x. This strategy reducesthe number of operations even more than usingHorner’s rule. Table 1 and 2 show the coefficientof the transformed polynomials for cubic and quin-tic order interpolation.4 PERFORMANCEANALYSISTable 3 shows the number of operations for thepolynomial evaluation for the classical 1-D cubicconvolution. Table 4 shows the results for ar-bitrary n-th order interpolation. As traditionalpiecewise n-th order convolution is separable (cas-caded convolution), the evaluation of N-D convo-lution takes N× 1-D convolution time.For 1-D cubic convolution the new method onlyneeds 16 instead of 26 operations, which means,that the proposed method manages to calculatethe same values with only 61.5% of operations.When assessing the whole interpolation task, wemust consider that the evaluation of the interpo-lation weights is one task, while the second task isto convolve the N-D function with these weights.The complexity of this convolution is given by(n + 1)N additions, multiplications and memoryaccesses. This basically means, that the convolu-tion dominates over the polynomial evaluation inhigher dimensions than 2.Therefore the achieved speedup, by transformingthe polynomials to [0, 1], for the complete interpo-lation task is higher for lower dimension interpo-lation, namely 1-D and 2-D.4.1 Real world experimentsA program was implemented in C++ to analysethe behavior of the new polynomial evaluation in areal world scenario. All calculation are done on aPentium III, 1 GHz. Of course the kernel functionscould be precalculated using lookup tables. Tokeep the accuracy at a high level we choose 10000grid points per unit length according to [Ostun97].Note that memory access is quite slow, thereforethe lookup table approach is slower than directcalculation for the cubic case. To get an averagevalue we looped over the calculation for 100 mil-lion times, which leads to the following results in2-D interpolation:Evaluation of the polynomials:using lookup tables 18988 msclassical method 13790 msnew method 7962 msWhole interpolation task:using lookup tables 51194 msclassical method 45996 msnew method 40168 msThe new polynomial evaluation only takes 57.7%of computation time in comparison to the clas-sical version, which matches quite well with thetheoretical value of 61.5%. The whole 2-D cubicconvolution task needs 87.3% in comparison to theclassical version.5 CONCLUSIONThis paper presents a novel method for efficientimplementation of cubic and higher order interpo-lation. The proposed polynomial evaluation needsonly about 60% of computation time for cubic andabout 65% for higher order polynomials in com-parison to the standard method. The whole imageinterpolation process is speeded up and takes onlyabout 87% for cubic convolution in comparison tothe standard method.6 ACKNOWLEDGEMENTSThanks to Christopher Zach for mathematicalsupport.A PROOFSLemma 1: The lowest order coefficients ai0 are0 for all polynomials fi(x), except for the firstpolynomial f0(x) where the coefficient is 1.Proof: The first polynomial f0(0) = 1⇒ a00 = 1.For all other polynomials fi(0) = 0 thereforeai0 = 0. Lemma 2: The highest order coefficients ainof pairwise symmetrical polynomials, fi(x) andfi+1(x), i odd, have the same values with oppositesign.Proof: The symmetrical polynomialscalled f(x) := fi(x) =∑nj=0 ajxj andg(x) := fi+1(x) =∑nj=0 bjxj have the prop-erty that f(x) = g(1− x). This leads to:f(x) = anxnn∑j=1ajxj = bn(1− x)nn∑j=1bj(1− x)jBy using factor comparison of an andbn, the coefficient for xn is given byanxn = bn(−x)n = (−1)nbnxn. Because nis odd (−1)n = −1 and this yields an = −bn. Lemma 3: All coefficients of order j ≤ n − 2of pairwise non symmetrical polynomials havesame the coefficients of odd order and the samecoefficients with opposite sign for even order.Proof: fi and fi+1, i ≥ 1 and i odd: Becauseof the continuous and smooth interpolationconstraints f (l)i (0) = (−1)lf (l)i+1(0), which yieldsail = (−1)lai+1l for l ≤ n− 2. Lemma 4: All coefficients anj with j ≤ n − 2 ofthe polynomial fn(x) are 0.Proof: Because of the continuous and smoothinterpolation constraints the last polynomialfn(x) must be continuous to the x-axis, thereforef(l)n (0) = 0 for l ≤ n − 2 which directly yieldsanj = 0 for j =≤ n− 2. Lemma 5: All coefficients a0j with odd or-der j with j ≤ n − 2 are 0 for the polynomialf0(x).Proof: f0(x) is defined as (equation 1)f0(x) =∑nj=0 a0j |x|j and is l-times continuousat x = 0. So limx→0− f (l)(x) = limx→0+ f (l)(x),which gives limx→0− f (l)(0) = a0l!=limx→0+ f (l)(0) = (−1)la0l. Therefore forall j ≤ n− 2 and j odd: a0j = 0. References[Keys81] R. G. Keys. Cubic convolution interpo-lation for digital image processing. IEEEai3 ai2 ai1 ai0f0 3/2 −5/2 0 1f1 −3/2 2 1/2 0f2 −1/2 1 1/2 0f3 1/2 −1/2 0 0Table 1: Coefficients of the transformed cubic convolution kernelsai5 ai4 ai3 ai2 ai1 ai0f0 −27/32 63/32 0 −17/8 0 1f1 27/32 −9/4 9/16 5/4 19/32 0f2 13/64 −19/64 −9/16 5/4 −19/32 0f3 −13/64 23/32 −9/32 −3/16 −3/64 0f4 3/64 −3/16 9/32 −3/16 3/64 0f5 −3/64 3/64 0 0 0 0Table 2: Coefficients of the transformed quintic convolution kernelsclassical method new method∗ + ∗ +calculation of distances 0 4 0 1polynomial evaluation 12 10 8 7overall 12 14 8 8Table 3: Number of needed operations for cubic convolution polynomial evaluationclassical method new method∗ + ∗ +calculation of distances 0 n+ 1 0 1polynomial evaluation n(n+ 1) n(n+ 1) n(n+ 3)/2 n2 − 3/2n+ 5/2overall n(n+ 1) (n+ 1)2 n(n+ 3)/2 n2 − 3/2n+ 7/2Table 4: Maximal number of needed operations for n-th order convolution polynomial evaluationclassical method new methodorder ∗ + ∗ + factor3 12 16 9 8 0.60715 30 36 20 21 0.62127 56 64 35 42 0.64179 90 100 54 71 0.6579Table 5: Maximal number of needed operations for 3rd − 9th order convolution polynomial evaluation.It could be possible that some coefficients are 0 or 1 by chance which decreases the number of neededoperations.Transactions on Acoustics, Speech, and Sig-nal Processing, 29(6):1153–1160, 1981.[Lehma99] T.M. Lehmann, C. Go¨nner, andK. Spitzer. Survey: interpolation methodsin medical image processing. IEEE Trans-actions on Medical Imaging, 18(11):1049 –1075, November 1999.[Meije99] E. Meijering, K. Zuidervel, andM. Viergever. Image reconstruction withsymmetrical piecewise nth-order polyno-mial kernels. IEEE Transactions on ImageProcessing, 8(2):192–201, February 1999.[Meije01] E. Meijering, W. Niessen, andM. Viergever. Quantitative evaluationof convolution-based methods for medi-cal image interpolation. Medical ImageAnalysis, 5:111–126, 2001.[Ostun97] J. L. Ostuni, A. K. S. Santha, V. S.Mattay, D. R. Weinberger, R. L. Levin, andJ. A. Frank. Analysis of interpolation effectsin the reslicing of functional MR images.J. Comp. Assisted Tomogr., 21(5):803810,1997.[Rifma73] S. S. Rifman. Digital rectification ofETRS multispectral imagery. Proc. Symp.Significant Results Obtained form the EarthResources Technology Satellite-1, Section B,NASA SP-327(1):1131–1142, 1973.[Unser99] Michael Unser. Splines: A perfect fit forsignal and image processing. IEEE SignalProcessing Magazine, 16(6):22–38, Novem-ber 1999.

English

Efficient implementation of higher order image interpolation

University of West Bohemia Digital Library

http://dspace5.zcu.cz/bitstream/11025/6248/1/J47.pdf

Efficient implementation of higher order image interpolation

Abstract

Similar works

Full text

Available Versions

DSpace at University of West Bohemia

University of West Bohemia Digital Library