2892

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 23, NO. 7, JULY 2014

A New Formula for Bivariate Hermite Interpolation on Variable Step Grids and Its Application to Image Interpolation Konstantinos K. Delibasis and Aristides Kechriniotis

Abstract— In this paper, we present a novel formula of the bivariate Hermite interpolating (BHI) polynomial in the case of support points arranged on a grid with variable step. This expression is applicable when interpolation of a bivariate function is required, given its value and the values of its partial derivatives of arbitrarily high order, at the support points. The proposed formula is a generalization of an existing formula for the bivariate Hermite polynomial. It is also algebraically much simpler, thus can be computed more efficiently. In order to apply Hermite interpolation to image interpolation, we simplify the proposed (BHI) to handle support points on a regular unit-step grid. The values of image partial derivatives are arithmetically approximated using compact finite differences. The proposed method is being assessed in a number of image interpolation experiments that include a synthetic image, for which the values of the partial derivatives are computed analytically, as well as a collection of images from different medical modalities. The proposed BHI with up to second-order image partial derivatives, outperforms the convolution-based interpolation methods, as well as generalized interpolation methods with the same number of support points that was compared with, in the majority of image interpolation experiments. The computational load of the proposed BHI is calculated and its behaviour with respect to its controlling parameters is investigated. Index Terms— Image interpolation, bivariate Hermite interpolation on variable step grids.

I. I NTRODUCTION

I

NTERPOLATION is required for every kind of geometric transformation of images, from simple rotation and resizing to complex elastic deformations. It is always desirable to minimize the artifacts of these geometric operations. Therefore, image interpolation is a subject that has attracted wide attention in literature, [1]–[3]. Reconstruction of medical images using the radon transform [1], as well as registration (spatial alignment) of two images of the same or different modality, requires image interpolation [4]. Interpolation-induced gray Manuscript received September 24, 2013; revised January 31, 2014 and April 16, 2014; accepted April 26, 2014. Date of publication May 7, 2014; date of current version May 27, 2014. Portions of the research in this paper use the CASIA-IrisV1 collected by the Chinese Academy of Sciences’ Institute of Automation (CASIA). The associate editor coordinating the review of this manuscript and approving it for publication was Prof. David Frakes. K. K. Delibasis is with the Department of Computer Science and Biomedical Informatics, University of Thessaly, Lamia 35100, Greece (e-mail: [email protected]). A. Kechriniotis is with the Department of Electronics, Technological Educational Institute of Lamia, Lamia 35100, Greece (e-mail: [email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TIP.2014.2322441

level errors of the transformed image affect the registration error in a number of imaging modalities [2]. A number of interpolation techniques have been proposed in the literature [1]. In this work we propose the use of bivariate Hermite interpolation (BHI) for image interpolation. The Hermite interpolation requires the values of the derivatives of the interpolated function, up to an arbitrary order, at each support point, in addition to the values of the interpolated function. In [7] and [8] the bivariate Hermite interpolation is discussed. In [9] the use of Hermite bivariate interpolation of the first order using 2 × 2 support points is proposed, whereas in [10] its use for surface reconstruction from shading is demonstrated. In [5] a formula is provided for the univariate Hermite interpolating polynomial, expressed as sums of derivatives of rational functions. In [6] we presented a much simpler algebraic new formula for the univariate Hermite interpolating polynomial that is expressed as sum of polynomials. However, the bivariate Hermite polynomial is not separable [32], [33], thus cannot be derived by applying the univariate one along the two dimensions. In [11], [12], formulas for the bivariate Hermite interpolating polynomial are provided, in the case of support points arranged on a non-equidistant grid. They are applicable in special conditions concerning the number of support points in each axis and the maximum order of partial derivatives for each support point. These formulas are algebraically very complex since they require the derivatives of rational functions. In this work, we propose a generalization of the formulas provided in [11], [12], removing the constraint of the aforementioned special conditions. Our proposed formula is also algebraically considerably simpler, since it constructs the Hermite polynomial as a sum of polynomials, rather than a sum of derivatives of rational functions as in [11], [12]. In addition, we make simplifying assumptions that allow the application of Bivariate Hermite Interpolation (BHI) to two-dimensional (2D) images. We combine the algebraic simplicity of the propose Hermite bivariate polynomial with a computational approach based on look up tables. In this way we achieve acceleration of the execution of the proposed BHI with selected settings, to a level that is comparable to convolution-based interpolation with separable kernels. The partial image derivatives that are required by the BHI are approximated using compact finite differences [13], rather than central finite differences [14]. This approach allows the use of

1057-7149 © 2014 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

DELIBASIS AND KECHRINIOTIS: NEW FORMULA FOR BHI ON VARIABLE STEP GRIDS

limited number of effective support points, which is a very desirable characteristic for any interpolation method. II. T HEORY OF H ERMITE I NTERPOLATION This section deals with the general case of BHI with support points on variable step grids. Image interpolation is a special case of interpolation with support points on unit step grid and it will be explored in the next section. A. Existing Theorems of Bivariate Hermite Interpolation on Grids In [11] the following bivariate generalization of Hermite’s interpolation was proven: Theorem 1: Let m be a positive integer that denotes order (i ,i ) of differentiation, p j11, j22 , 0 ≤ i 1 , i 2 < m, 1 ≤ j1, j2 ≤ n be given real numbers; let (λ j1 , μ j2 ) be n 2 distinct points for 1 ≤ j1 , j2 ≤ n. Then there exists a unique polynomial p (x 1 , x 2 ) with degx1 p < nm, degx2 p < nm such that (i1 ,i2 ) ∂ i1 +i2 i1 i2 ( p(λ j1 , μ j2 )) = p j1 , j2 , 0 ≤ i 1 , i 2 < m, 1 ≤ j1 , ∂ x1 ∂ x2 j2 ≤ n.

This polynomial is given by:

p (x 1 , x 2 ) =

[L 1 (x 1 ) L 2 (x 2 )]m

×

[(m − 1)!]2 n m−1 n    m−1  j1 =1 j2 =1 i1 =0 i2 =0

(i ,i )

(i ,i )

M j11, j22 (x 1 , x 2 ) p j11, j22 ,

(1)

where L 1 (x 1 ) =

n  



x 1 − λ j1 , L 2 (x 2 ) =

j1 =1

n  

 x 2 − μ j2 ,

j2 =1

     1 m −1 m −1 d m−1−i1 (i1 ,i2 ) M j1 , j2 (x 1 , x 2 ) = i2 i1 ds m−1−i1 A j1 (s) s=λ j 1  m−1−i2   d 1 · dt m−1−i2 B j2 (t) t =μ j 2

with

m L 1 (s)  , A j1 (s) = (s − x 1 )  s − λ j1

m L 2 (t)  . B j2 (t) = (t − x 2 )  t − μ j2

In [12], Theorem 1 is generalized so that the maximum order of the partial derivative with respect to x 1 at any point λ j1 , μ j2 depends only on j1 and the maximum order of the partial derivative with respect to x 2 depends only on j2. Using [12] we conclude that the error of the interpolation formula of Theorem 1 is the following: R (x 1 , x 2 ) = f (x 1 , x 2 ) − p (x 1 , x 2 ) [L 1 (x 1 )]m ∂ nm [L 2 (x 2 )]m ∂ nm = f (x 1 , ξ2 ) nm ‘ f (ξ1 , x 2 ) + (nm)! ∂ x 1 (nm)! ∂ x 2nm [L 1 (x 1 ) L 2 (x 2 )]m ∂ nm+nm f (ξ1 , ξ2 ) , − (2) ∂ x 1nm ‘ ∂ x 2nm [(nm)!]2

2893

   for a point  (ξ1 , ξ2 ) ∈ min x 1 , λ j1 , . . . , λ jn , max x 1 , λj1 , . . . , λ jn × min x 2 , μ j1 , . . . , μ jn , max x 2 , μ j1 , . . . , μ jn . An integral formula of the error R(x 1 , x 2 ) is given in [11]–[12]. B. The Proposed Bivariate Hermite Interpolation In this subsection we derive a new formula for the bivariate Hermite polynomial that: a) It is algebraically simpler than (1), since it does not contain derivatives of rational functions, thus it is calculated more efficiently and b) It can be applied to more general circumstances, since the maximum order of the partial   derivatives at any point λ j1 , μ j2 depends on both j1 and j2 . Before presenting the Theorem for the proposed expression of the BHI polynomial, we need the following definitions.   Let us suppose that for each support position λ j1 , μ j2 , 1 ≤ j1 ≤ n 1 , 1 ≤ j2 ≤ n 2 two non-negative integers k j1, j2 , m j1 , j2 are given, such that the value of all partial   derivatives  of a bivariate function f (x 1 , x 2 ), ∂ i1 +i2 f λ , 1 ≤ i 1 < k j1 , j2 , 1 ≤ i 2 < m j1 , j2 are , μ j j i1 i2 1 2 ∂ x1 ∂ x2

prescribed. The value of the aforementioned partial derivatives (i ,i ) of f will be denoted as p j11, j22 , (the two subscripts denote the index of the support position and the two superscripts indicate the order of partial derivative with respect to the first and second coordinate). The proposed BHI, constructs a unique polynomial p (x 1 , x 2 ), which interpolates the value of the function f (x 1 , x 2 ) as well as the given values of its partial derivatives of all given orders, at each support position. We can formalize the proposed BHI as following. Let us define the following polynomials    x 1 − λh 1 k j1 , j2 K j1 , j2 (x 1 ) =  k j , j , 1 2 λ − λ j h h 1  = j1 1 1  m j , j  x 2 − μh 2 1 2 L j1 , j2 (x 2 ) = . (3)  m μ j 2 − μh 2 j 1 , j 2 h = j 2

2

Let us also define the rows with dimensions 1 × k j1 , j2 and 1 × m j1 , j2 respectively:   k j , j −1 1 2 x −λ x −λ X j1 , j2 (x 1 ) = K j1 , j2 (x 1 ) 1, 1 1! j1 , . . . , 1 j1  , −1 ! k j1 , j2

 Y j1 , j2 (x 2 ) = L j1 , j2 (x 2 ) 1,

x 2 −μ j2 1! ,



...,

m j

−1 1 , j2

x 2 −μ j2   m j1 , j2 −1 !

(4a) . (4b)

Let us denote by V the vector space spanned by the polyno    x −λ

i1

x −μ

i2

2 j2 K j1 , j2 (x 1 ) L j1 , j2 (x 2 ) , 1 ≤ j1 ≤ n 1 , mials 1 i1 !j1 i2 ! 1 ≤ j2 ≤ n 2 , 0 ≤ i 1 < k j1 , j2 , 0 ≤ i 2 < mj1 , j2 . Finally, let us define for each support position λ j1 , μ j2 the k j1 , j2 × k j1, j2 and m j1 , j2 × m j1 , j2 lower unit triangular matrices  j1 , j2 and  j1 , j2 , as (5) and (6),  shown at the top of the next page. Theorem 2: Let λ j1 , μ j2 , 1 ≤ j1 ≤ n 1 , 1 ≤ j2 ≤ n 2 be n 1 × n 2 discrete support positions that lie on a variable step grid with distinct λ j1 and distinct μ j2 coordinates. (i ,i ) If k j1 , j2 , m j1 , j2 are positive integers and p j11, j22 , 0 ≤ i 1
0). Non-separability does increase

TABLE II C OMPARISON OF THE P ROPOSED BHI U NDER 3 D IFFERENT S ETTINGS W ITH S TATE OF THE A RT M ETHODS , IN T ERMS OF RMSE AND NSD FOR THE S YNTHETIC I MAGE I NTERPOLATION IN THE C ASE OF 4 AND 9 ROTATIONS

significantly the computational work load in the case where the fractional parts of the required image coordinates λreq and μreq are constant for any pixel location (e.g. image upsampling). Such an example is provided in [21], where it is stated that in the case of n point kernel applied in two dimensions, 2n evaluations are needed, compared to n 2 evaluations for a non-separable kernel. However more elaboration is needed for a more general geometric transformation, like the affine one, in which the x 1 and x 2 depend on the pixel position. In this case, in order to exploit the separability property of a kernel, the transformation needs to be decomposed into a series of one-dimensional transformations [22], (a relevant example of rotation decomposition into three one-dimensional translation or skew transformations is given in Gotchev et al. [23, page 321]). It has been shown that any 2D affine transformation can be decomposed into 3 one-dimensional transformations [22], [24]. Therefore, an n point separable kernel has to be applied 3 times in the case of a 2D affine geometric transformation. Equivalently, no computational gain is achieved for n < 4, (n the number of nominal support points per image axis). In this work, the proposed BHI method is employed for n = 2 (2 × 2 nominal support points), since as it will be shown in the Results section, in most of the experiments it outperforms the state of the art methods in comparison with same number of effective support points, including to O-MOMS and B-spline of 7th degree. The calculation of discrete derivatives using compact differences possesses the separability property, [17], [18]. Furthermore, if the 2D geometric transformation is not separable into a series of one-dimensional transformations (such as a non-linear or elastic transformation), a separable kernel offers no computational advantage.

DELIBASIS AND KECHRINIOTIS: NEW FORMULA FOR BHI ON VARIABLE STEP GRIDS

2897

Fig. 3. The synthetic image, one mammographic and one iris image used in the interpolation experiments with their portions used in Fig. 4 and 5 outlined.

It would seem convenient to express the proposed BHI as a two-dimensional kernel. For instance, in the case of the univariate Hermite interpolating polynomial, described in [5] and more recently in [6], the corresponding kernels have been produced for different number of support points [26], approximating discrete derivatives by central differences. In the current work however, the non-separable bivariate Hermite interpolating polynomial causes the corresponding Hermite kernels to be two-dimensional, thus they would have a large number of branches. In addition, the use of compact finite differences necessitates one such kernel for each image partial derivative, thus the kernel option becomes less practical. A single kernel could be generated using central differences instead of compact differences for derivative approximation, however the performance of the BHI would be inferior and the number of the required effective support points would increase. Therefore, the direct implementation of (8) is preferred. IV. R ESULTS A. Data The proposed evaluation scheme has been applied to one synthetic image, as well as to a number of images from different medical modalities. More specifically, the synthetic image with dimensions 512 × 512 pixels was constructed as a radial sinusoidal chirp with higher frequency at the image center (see Fig. 3), similar to the one used in the work of Thevenaz et al [25]. The synthetic image was selected since its partial derivatives can be calculated analytically, as well as approximated arithmetically, thus allowing to study the behavior of the proposed BHI in more detail. From the medical modalities, the following were used: 8 images 1024 × 1024 from the Mini-MIAS mammography database [27] and 8 images from the Iris database v1 with dimensions 320 × 200 [28]. One iris and one mammographic image are shown in Fig. 3, each with a portion indicated to be used for error visualization in Figs. 4 and 5. Furthermore, 10 transverse slices at the level of lower abdomen from a Computer tomography study, denoted as ‘CT’, with dimensions of 512 × 512 and 10 sagittal slices, 256 × 256 from an MRI brain study (‘MRI’). B. Scheme for the Evaluation of BHI The BHI is quantitatively evaluated using a series of affine geometric transformations that are applied sequentially on any given image I0 and result in a final image I f that ideally

Fig. 4. The absolute difference between the ideal image and the result of 9 rotations of the synthetic image, for the proposed BHI using n = 2, M = 2, with analytic derivatives, as well as with arithmetic derivative approximation (d1 = 10, d2 = 6). Other state of the art interpolation methods are also included.

should be identical to the original image. In this work, we consider the following: A) Series of 2, 4 and 8 pairs of random forward – inverse affine transformation: a random forward affine transformation is applied to a given image I0 to produce the image I1 , followed by the application of its inverse applied to the resulting I1 . A series of pairs of forward – inverse transformations are applied successively (see [2]–[29]). Each affine transformation consists of image rotation and scaling with respect to its centre. The angle and scaling factors are selected randomly in the range [0, π/3] and [0.9, 1.1]. B) Series of 2, 4 and 9 randomly selected successive rotations, which add to a total of 2π are applied to a given image to produce the final image I f , in an approach similar to [1], [2], [25, subsection III-A], [37 subsection 3.4]. The geometric transformations are implemented using the proposed BHI as well as the other interpolation methods under comparison. The root mean square error (RMSE) between I f and I0 was used as an objective measure of the quality of the investigated interpolation method [2]. We have also calculated the number of sites of disagreement (NSD) (see [30], [31]) between I f and I0 , defined as the number

2898

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 23, NO. 7, JULY 2014

Fig. 5. The resulting absolute difference between the ideal and the interpolated portion of an iris image (right) and a mammography image (left), for the proposed BHI and other state of the art interpolation methods, for 8 pairs of forward-inverse affine transforms. Black indicates zero difference. The interpolation methods are arranged in 2 × 4 as shown in the central part of the Figure. The image portions are shown on the original images in Fig. 3.

of pixels with values that differ from the correct ones by more than a predefined threshold. This threshold was set to 10 for the image interpolation experiments in this work. In order to avoid border effects, the images were mirrored in both dimensions ([2, subsection 4.1], [32, subsection 2.1]). C. Other Interpolation Techniques in Comparison A number of different interpolation techniques were selected to compare against the proposed BHI interpolation, including linear interpolation and Lagrange interpolation [33] for 4 and 6 support positions, implemented as in [1]. From the family of cubic interpolation kernels with 4 support positions, defined in [1,(23)], we used the ones generated by setting the kernel free parameter a to −1/2 which corresponds to Catmull-Rom interpolation [34], −3/4 and −1. We also included comparison with the Mitchell-Netravali (MN) kernel [35] (as defined in Eq.(26) from [1] using b = c = 1/3), the cubic polynomial interpolation using 6 and 8 support positions using the kernel, provided in [1 (24)] and the 4-point Lanczos interpolation [15 (16.61)]. From the methods of generalized interpolation, which are considered state of the art, we included B-spline interpolation ([20], [36]), as well as the optimal smallest-support support functions (O-MOMS) [37] of up to the 7th degree. This degree was selected since the kernels use the same number of support points per image dimension (8 × 8) with the effective support

points of the proposed BHI for n 1 = n 2 = 2, d1 = 10, d2 = 6, 8. The BHI for n 1 = n 2 = 2, with d1 ≤ 8, d2 ≤ 6, involves effective support points in a 6 × 6 neighborhood, thus it may be compared to the B-spline and O-MOMS of the 5th degree. Finally, we included the FT-based two stage resampling with magnification factor of 2, proposed in [32], in conjunction with the cubic polynomial 6- and 8-point interpolation, as well as in conjunction with the B-splines and the O-MOMS. All interpolation methods in comparison were implemented using Matlab. D. Comparative Results 1) Synthetic Image: Table II shows the interpolation error in terms of RMSE and NSD for 4 and 9 successive rotations of the synthetic image. The time required for one rotation using the aforementioned Matlab implementations is also given. To facilitate comparison, the interpolation methods are grouped according to their number of support points (see [2, p119]). Three entries for the proposed BHI are shown, using n 1 = n 2 = 2, M = 2. In the first BHI entry, the image partial derivatives are computed from the analytic expression of the image data, thus the arithmetic approximation is not applied. Therefore, only image information on the 2 × 2 nominal support points is used in the interpolation process. The other two BHI entries are using the arithmetic image derivative approximation with d1 = 8, d2 = 6 (thus using 6 × 6 effective support points)

DELIBASIS AND KECHRINIOTIS: NEW FORMULA FOR BHI ON VARIABLE STEP GRIDS

2899

TABLE III T HE RMSE OF THE P ROPOSED BHI AND THE O THER I NTERPOLATION M ETHODS , AVERAGED OVER THE AVAILABLE I MAGES OF E ACH M ODALITY FOR THE

S ERIES OF F ORWARD – I NVERSE A FFINE T RANSFORMATION E XPERIMENTS , G ROUPED BY THE N UMBER OF S UPPORT P OINTS . T HE P ROPOSED BHI W ITH 6 × 6 E FFECTIVE S UPPORT P OINTS HAS PARAMETERS n 1 = n 2 = 2, M = 2, d1 = 8, d2 = 6, W HEREAS FOR 8 × 8 E FFECTIVE S UPPORT P OINTS BHI HAS PARAMETERS n 1 = n 2 = 2, M = 2, d1 = 10, d2 = 6

and for d1 = 10, d2 = 6 (8 × 8 effective support points). The entries referring to the two-stage resampling method [32] are indicated by “FT”. The proposed BHI with analytically calculated image derivatives, although it uses only 2 × 2 support points, achieves the smallest RMSE and NSD of all the interpolation methods in comparison, including the two implementations of the proposed BHI that use arithmetic image derivative approximations. The BHI with 6 × 6 and 8 × 8 effective support points outperform the state of the art methods in comparison with the same number of support points (including the B-spline and the O-MOMS of the 5th and 7th degree respectively). It can also be observed that the BHI with 6 × 6 effective support points outperforms all the interpolation methods, including those with 8 × 8 support points as well, including the methods combined with the FT-based two stage resembling. The FT resampling causes a marked improvement in the case of cubic polynomial interpolation, or 3rd degree B-spline and more modest improvement when combined with the 7th degree O-MOMS. Furthermore, it is verified that the generalized interpolation methods (the B-spline and the O-MOMS of various degrees) clearly outperform the classic interpolation methods (e.g. the 3rd degree B-splines with 4 × 4 support points, achieves smaller RMSE than the cubic polynomial interpolator with 8 × 8 support points). Fig. 4 visualizes the resulting interpolated image after 9 rotations of the synthetic image for selected interpolation techniques, including BHI with analytic and numeric derivative calculation and the state of the art B-spline and O-MOMS

interpolators. The images appear to correlate well with the numerical results in Table II. 2) Real Images: Tables III and IV show the RMSE and the NSD of the interpolation methods, averaged over the available images, from mammography, iris, CT and MRI studies, for 2, 4 and 8 pairs of forward – inverse affine transformations. Table V shows the RMSE for 2, 4 and 9 successive rotations. The proposed BHI for n 1 = n 2 = 2, M = 2 is evaluated for two different settings of the order of accuracy of the first and second partial derivatives: a) d1 = 8, d2 = 6 that corresponds to 6 × 6 effective support points and b) d1 = 10, d2 = 6 that corresponds to 8 × 8 effective support points. To facilitate comparison, the interpolation methods are grouped according to their number of support points [2, p119]. It can be observed that in 21 out of 24 interpolation experiments in Tables III and V, the proposed BHI with d1 = 8, d2 = 6 outperforms in terms of RMSE the state of the art interpolation with the same number of support points (6 × 6) or less (including the B-spline and the O-MOMS of the 5th degree and the cubic interpolating polynomial, even combined with the FT-based two-stage resampling). In terms of NSD, the proposed BHI with d1 = 8, d2 = 6 outperforms the other interpolation methods with same or less support points (6 × 6) in 8 out of the 12 experiments of Table IV. Furthermore, in 19 out of the 24 interpolation experiments reported in Tables III and IV, the proposed BHI with d1 = 10, d2 = 6 achieves lower RMSE than any of the state of the art interpolation methods with the same number of support points (8 × 8) or less. In terms of NSD, the proposed BHI with d1 = 10, d2 = 6 outperforms the other interpolation methods with same or

2900

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 23, NO. 7, JULY 2014

TABLE IV T HE N UMBER OF S ITES OF D ISAGREEMENT (NSD) FOR THE I NTERPOLATION M ETHODS AND THE E XPERIMENTS OF TABLE III

TABLE V T HE RMSE FOR THE ROTATION E XPERIMENTS ( SEE C APTION OF TABLE III)

less support points (8 × 8) in 9 out of the 12 experiments of Table IV. As expected, the generalized convolution methods (O-MOMS and B-spline) consistently outperform the traditional convolution methods (the Lanczos, MN, Lagrange and cubic polynomial interpolator with various settings). The FT-based two-stage resampling improves significantly the performance of the 6-point and 8-point cubic polynomial interpolation, both in terms of RMSE and NSD. Significant

improvement is also achieved when the FT-based two-stage is combined with B-spline of 3rd degree, whereas marginal improvement is observed with the combination of the FT two-stage with the 7th degree O-MOMS. In order to assess the performance of the proposed BHI, we present in Fig. 5 the resulting difference images for 8 affine transformation pairs for a portion of an iris image and for a portion of one of the mammography images (image portions are shown in Fig. 3). The BHI with 8 × 8 support points and 7 state of

DELIBASIS AND KECHRINIOTIS: NEW FORMULA FOR BHI ON VARIABLE STEP GRIDS

2901

Fig. 7. The RMSE of the BHI for the synthetic image with 4 rotations, as function of n and M. The RMSE of the 7th degree O-MOMS is also shown. Fig. 6. The achieved RMSE of various state of the art interpolation methods grouped by the number of support points, for the 4-rotation experiment of the Mammography images.

the art interpolation methods are arranged in 2 columns and 4 rows, as shown in the central part of Fig. 5. The visual appearance of the difference images correlates well with the numerical results in Table III. Fig. 6 summarizes the RMSE for various state of the art interpolation methods, grouped by the number of their support points (the effective number of support points was used for the BHI). The 4-rotation experiment using the Mammography images was selected, as a typical example of the majority of experiments in which the proposed BHI using 8 × 8 support points outperforms the other methods. V. D ISCUSSION A. Computational Complexity of the Proposed BHI for Image Interpolation Considering the computational implementation of BHI using look up tables, as described previously, the calculations of (7) involve one multiplication of the 1 × (M + 1) row, times the (M + 1) × (M + 1) matrix, times one column (M +1)×1, which add up to (2M +1)(M +1)+ M +1+ M = (2M + 1)(M + 2) arithmetic operations for each support point. Thus, the complexity for n×n nominal support points becomes n 2 (2M + 1)(M + 2). The approximation of the image partial derivatives imposes a computational cost, which depends on the order of accuracy d1 , d2 . Setting d1 equal 6, 8 and 10, the size of the row convolution kernel at the right hand side of Eq.(12) is 5, 5 and 7 respectively. The number of arithmetic operations introduced by the left hand side of (10) and (11) is O(4N), O(8N) and O(8N) respectively, as described in detail in [19], [20]. Concerning the generalized interpolation (B-splines and O-MOMS), the cost of IIR prefiltering for each image pixel is obtained as described in detail in [20]. The cost of the subsequent convolution per image pixel with an n×n separable kernel is 3(2n − 1) arithmetic operations in the case of an affine geometric transformation that can be decomposed into 3 one-dimensional transformations (see par. 2.2.3), or 2n 2 − n operations in the case of a non-decomposable geometric

transformation, in which case the convolution step needs to be performed in two dimensions. The proposed BHI can be computationally more demanding than the generalized interpolation methods, for separable affine geometric transformations, but it presents the following advantages: a) it achieves lower interpolation RMSE, as well as NSD in the majority of the interpolation experiments, compared to kernels with the same number of support points and b) it is the method of choice when the values of the derivatives are not arithmetically approximated from the image data, but originally acquired (analytically or by measurements), as it has been shown in the results with the synthetic image. For geometric transformations that cannot be decomposed in series of 1D transformations, the difference in computational complexity becomes less significant. B. The Behaviour of BHI With Respect to Its Controlling Parameters As it has been mentioned previously, the controlling parameters of the BHI are the nominal number of support points per image axis n(n × n points in two dimensions), the maximum order of partial differentiation M and the order of accuracy of the 1st and 2nd image derivatives d1 , d2 respectively, which define the number of effective support points. We can rewrite the error expression (2) for Theorem 1 for any point (ξ1 , ξ2 ) ∈ [ j1, j1 + n − 1] × [ j2, j2 + n − 1], under the assumptions presented in section 3: f (x 1 , x 2 ) − p (x 1 , x 2 ) (x 1 − j1) M+1 K j1 (x 1 ) ∂ n(M+1) = f (ξ1 , x 2 ) n(M+1) (n (M + 1))! ∂x 1

(x 2 − j2 ) M+1 L j2 (x 2 ) ∂ nm + f (x 1 , ξ2 ) n(M+1) (n (M + 1))! ∂x 2

−

((x 1 − j1) (x 2 − j2 )) M+1 K j1 (x 1 ) L j2 (x 2 )

×

[(n (M + 1))!]2 ∂ 2n(M+1) n(M+1)‘

∂ x1

n(M+1)

∂ x2

f (ξ1 , ξ2 ) ,

(12)

Fig. 7 investigates the behaviour of the proposed BHI with respect to its most important controlling parameters: n, M, whereas d1 and d2 were set equal to 10 and 6 respectively

2902

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 23, NO. 7, JULY 2014

d1 ∈{6,8,10}, d2 ∈{ 4,6,8} . The rest of the BHI parameters were set to n = 2 (2 × 2 nominal support points) and M = 2, as discussed before. The RMSE of O-MOMS of 7th degree (8 × 8) is also included. It can be observer that a significant RMSE improvement occurs when increasing the order of accuracy of the 1st image derivative d1 from 6 to 8 and a more marginal improvement by setting d1 = 10. This is expected, since increasing d1 results in increasing the number of effective support points. Similar behaviour is observed in the majority of the interpolation experiments. Increasing d2 from the value of 4 up to 8 also causes a reduction of the RMSE, when d1 = 8 or 10, although less noticeable. The RMSE of O-MOMS of 7th degree (one of the best performing methods in comparison), is also included in the Figure. Similar behaviour has been observed in most of the experiments with the other image modalities.

Fig. 8. The RMSE of the proposed BHI for an image of the mammography database and 1 pair of forward - inverse affine transformations, as a function of 1st and 2nd derivative accuracy order.

(d2 is not applicable for M = 1). The synthetic image under 4 rotations with arithmetic approximation of the partial derivatives was used. It can be observed that increasing n while keeping M = 1, results in significant error improvement for M = 1 (Fig. 7). This behavior is expected from the expression of the interpolation error (12), since all the tree error terms contain n! in the denominator. On the other hand, each error term also contains the value of the (n(M +1)) order of the function partial derivatives, which cannot be predicted, which may cause the behavior for M = 2. Furthermore, all BHI with M = 2 outperform the BHI with M = 1, at the expense of slight increment to its computational complexity. This behavior is also expected from (12). As M increases, the denominator of each error term increases in a factorial manner. If n = 2 the nominator of each fraction decreases exponentially. If n > 2 the nominator increases exponentially. Therefore, the overall value of the fraction of each error term decreases. In Fig. 8 the RMS error of the BHI is calculated for an image of the mammography database and 1 pair of forward - inverse affine transformations, using all combinations of

∂ i1 +i2



x 1 − λt1

s1 

x 2 − μt 2 s1 !s2 !

VI. C ONCLUSIONS A new formula for the bivariate Hermite interpolation for grids (BHI) for support points arranged on a non-uniformly spaced grid is derived in this work. The formula is more general and algebraically simpler than the previously proposed ones. A number of simplifying assumptions are made to achieve an efficient computational implementation for image interpolation. The required image partial derivatives are calculated using compact differences. The controlling parameters of the BHI are identified and their effect on the proposed method is studied, both theoretically and experimentally. The computational complexity of the proposed method is also studied in detail. The proposed BHI is compared against a number of state of the art interpolation techniques, using extensive experiments for a synthetic image with known partial derivatives, as well as a number of images from different medical imaging modalities. Results show that the proposed BHI outperforms the state of the art interpolation techniques that it was compared with, in terms of RMS error and NSD in the majority of the interpolation experiments, while using the same number of support points. 

s2

K j1 , j2 (x 1 ) L j1 , j2 (x 2 )   ∂ x 1i1 ∂ x 2i2 (x 1 ,x 2 )= λt1 ,μt2 ⎧ 0 i f (t1 , t2 ) = ( j1, j2 ) ⎪ ⎪ ⎪ ⎨0 i f (t , t ) = ( j , j ) and (i < s ori < s ) 1 2 1 2 1 1 2 2 =           ⎪ i1 i2 (i −s ) (i −s )  ⎪ ⎪ K j1 , j2 1 1 λ j1 L j1 , j2 2 2 μ j2 , i f (t1 , t2 ) = ( j1 , j2 ) , s1 ≤ i 1 , s2 ≤ i 2 ⎩ s1 s2 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

(0,0)

(0,1)

c j1 , j2

c j1 , j2

  k j1 , j2 −1,0 c j1 , j2

  k j1 , j2 −1,1 c j1 , j2

(1,0) c j1 , j2

(1,1) c j1 , j2

.. .

...

 0,m j

c j1 , j2 1



, j2 −1

⎤

⎡

(0,0)

p j1 , j2

(0,1)

p j1 , j2

...

 0,m j

p j1 , j2 1

, j2 −1



(A1)

⎤

⎥ ⎢  ⎥  ⎥ ⎢ ⎥ 1,m j1 , j2 −1 (1,1) ⎥  ⎥  −1 ⎢ p(1,0)  T p j1 , j2 p j1 , j2 j1 , j2 ⎥=  j , j ⎢ ⎥  j . j −1 1 2 1 2 ⎥ ⎢ ⎥ .. .. .. ⎥ ⎢ ⎥ . . . ⎣  ⎦    ⎦   kj1 , j2 −1,m j1 , j2 −1 k j1 , j2 −1,0 k j1 , j2 −1,1 k j1 , j2 −1,m j1 , j2 −1 p j1 , j2 . . . c j1 , j2 p j1 , j2 . . . p j1 , j2 (A2)   1,m j1 , j2 −1 c j1 , j2

DELIBASIS AND KECHRINIOTIS: NEW FORMULA FOR BHI ON VARIABLE STEP GRIDS

⎡ ⎢ ⎢ ⎢  j1 , j2 ⎢ ⎢ ⎢ ⎣ ⎡ ⎢ ⎢ ⎢ ⎢ =⎢ ⎢ ⎢ ⎢ ⎣

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

(0,0)

(0,1)

c j1 , j2

c j1 , j2

c j1 , j2

c j1 , j2 .. .

(1,0)

  k j1 , j2 −1,0 c j1 , j2

  k j1 , j2 −1,1 c j1 , j2

k j , j −1 1 2 k j , j −1 ∂ x1 1 2

.. .

k j , j −1 1 2 k j , j −1 ∂ x1 1 2

∂

⎡

⎢ ⎢ ⎢ =⎢ ⎢ ⎢ ⎣

...



 k j , j −1,m j1 , j2 −1 c j1 ,1j2 2

k ∂ j1 , j2 k j , j −1 1 ∂ x1 2 ∂ x2

  p λ j1 , μ j2   ∂2 ∂ x 1 ∂ x 2 p λ j1 , μ j2 k ∂ j1 , j2 k j , j −1 1 ∂ x1 2 ∂ x2

(0,1)

p j1 , j2

p j1 , j2

p j1 , j2 .. .

...

(1,1)

  k j , j −1,1 p j1,1j2 2

···

 0,m j

, j2 −1

  1,m j1 , j2 −1 ,

.. .

k j , j +m j , j −2 1 2 1 2 k j , j −1 m j , j −1 1 2 ∂ x1 ∂ x2 1 2

∂

 p λ j1 , μ j2

  p λ j1 , μ j2   p λ j1 , μ j2

m j , j −1 1 2 m j , j −1 ∂ x2 1 2 m ∂ j1 , j2 m j , j −1 ∂ x1 ∂ x2 1 2

···



p j1 , j2 1

.. .

∂

...

⎤

  p λ j1 , μ j2   p λ j1 , μ j2

m j , j −1 1 2 m j , j −1 ∂ x2 1 2 m ∂ j1 , j2 m j , j −1 ∂ x1 ∂ x2 1 2

∂

...

  p λ j1 , μ j2 . . .

p j1 j2 ...

⎥ ⎥ ⎥  ⎥ j ,j T 1 2 ⎥ ⎥ ⎦

  p λ j1 , μ j2 . . .

.. .

k j , j +m j , j −2 1 2 1 2 k j , j −1 m j , j −1 ∂ x1 1 2 ∂ x2 1 2

∂

 p λ j1 , μ j2

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

(A4)

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

⎤

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

 k j , j −1,m j1 , j2 −1 p j1,1j2 2

VII. A PPENDIX A Using the Leibniz’s rule for derivatives, we easily get the following Lemma: Lemma 1: For 1 ≤ j1 , t1 ≤ n 1 , 1 ≤ j2 , t2 ≤ n 2 , 0 ≤ i 1 , s1 < k j1 , j2 , 0 ≤ i 2 , s2 < m j1 , j2 as (A1), shown at the bottom of the previous page. Proof of Theorem 1: It can be observed that (8) can be equivalently written in the following form:  i  i , j2 −1 m j1 , j2 −1 n1  n 2 k j1  (i ,i ) x 1 −λ j1 1 x 2 −μ j2 2 p (x 1 , x 2 ) = c j11, j22 i1 ! i2 ! j1 =1 j2 =1 i1 =0

⎤



  p λ j1 , μ j2   ∂2 ∂ x 1 ∂ x 2 p λ j1 , μ j2

p j1 , j2

  k j , j −1,0 p j1,1j2 2

, j2 −1

∂ ∂ x2

  p λ j1 , μ j2

(1,0)

 1,m j

c j1 , j2 1 .. .

.. .

(0,0)

, j2 −1

∂ ∂ x2

  p λ j1 , μ j2

  p λ j1 , μ j2   ∂ ∂ x 1 p λ j1 , μ j2

 0,m j

c j1 , j2 1

(1,1)

  p λ j1 , μ j2   ∂ ∂ x 1 p λ j1 , μ j2

∂

...

2903

i2 =0

K j1 , j2 (x 1 ) L j1 , j2 (x 2 ), where as (A2), shown at the bottom of the previous page. Now by calculating the derivative of order (i 1, i 2 ) of the  polynomial p at λ j1 , μ j2 and using (A1) of Lemma 1 we obtain:    i2 i1     ∂ i1 +i2 i2 (k1 ,k2 ) i 1 = p λ , μ c j1 j2 j1 , j2 i1 i2 k k 1 2 ∂ xi ∂ x2 k1 =0 k2 =0 (i −k )    (i −k )    (A3) K j1 , j2 1 1 λ j1 L j1 , j2 2 2 μ j2 for all 1 ≤ j1 ≤ n 1 , 1 ≤ j2 ≤ n 2 , 0 ≤ i 1 < k j1, j2 ,0 ≤ i 2 < m j1 , j2 . The system of equations (A4) can be rewritten as (A4), shown at the top of the page. Substituting (A2) in (A4), we get the equation shown at the top of the page below eq.(A4)

Thus, we derive that

∂ i1 +i2 i i ∂ x i 1 ∂ x 22 ≤ j2

  p λ j1 , μ j2 =

(i ,i )

p j11, j22

for

all 1 ≤ j1 ≤ n 1 , 1 ≤ n 2 , 0 ≤ i 1 < k j1 , j2 , 0 ≤ i 2 < m j1 , j2 . Finally, for the uniqueness of the polynomial p it  is sufficient to prove that the polynomials    x 1 −λ j1 i1 !

i1

x 2 −μ j2 i2 !

i2

K j1 , j2 (x 1 ) L j1 , j2 (x 2 ) , for all i 1 , i 2 , j1 , j2 are linearly independent: We start from p (x 1 , x 2 ) =  i  i , j2 −1 , j2 −1 m j1 n 2 k j1 n1   (i ,i ) x −λ 1 x 2 −μ j2 2 c j11, j22 1 i1 !j1 K j1 , j2 (x 1 ) i2 ! j1 =1 j2 =1 i1 =0

i2 =0

L j1 , j2 (x 2 ) = 0. Then since the polynomial p (x 1 , x 2 ) is identically equal to zero, we have that (A2 ) states by (i ,i ) p j11, j22 = 0 for all 1 ≤ j1 ≤ n 1 , 1 ≤ j2 ≤ n 2 , 0 ≤ i 1 < k j1 , j2 , (i ,i )

0 ≤ i 2 < m j1 , j2 . Therefore from (A2 ) we get c j11, j22 = 0. ACKNOWLEDGMENT

Portions of the research in this paper use the CASIA-IrisV1 collected by the Chinese Academy of Sciences’ Institute of Automation (CASIA). R EFERENCES [1] T. M. Lehmann, “Survey: Interpolation methods in medical image processing,” IEEE Trans. Med. Imag., vol. 18, no. 11, pp. 1049–1075, Nov. 1999. [2] E. H. W. Meijering, W. J. Niessen, and M. A. Viergever, “Quantitative evaluation of convolution-based methods for medical image interpolation,” Med. Image Anal., vol. 5, no. 2, pp. 111–126, 2001.

2904

[3] E. H. W. Meijering, W. J. Niessen, J. P. W. Pluim, and M. A. Viergever, “Quantitative comparison of sinc-approximating kernels for medical image interpolation,” in Medical Image Computing and ComputerAssisted Intervention (Lecture Notes in Computer Science), vol. 1679, C. Taylor and A. Colchester, Eds. Berlin, Germany: Springer-Verlag, 1999, pp. 210–217. [4] J. Tsao, “Interpolation artifacts in multimodality image registration based on maximization of mutual information,” IEEE Trans. Med. Imag., vol. 22, no. 7, pp. 854–864, Jul. 2003. [5] I. S. Berezin and N. P. Zhidkov, Computing Method. New York, NY, USA: Pergamon, 1973, ch. 8, sec. 9. [6] K. Delibasis, A. Kechriniotis, and N. Assimakis, “New closed formula for the univariate hermite interpolating polynomial of total degree and its application in medical image slice interpolation,” IEEE Trans. Signal Process., vol. 60, no. 12, pp. 6294–6304, Dec. 2012. [7] R. A. Lorentz, “Multivariate Hermite interpolation by algebraic polynomials: A survey,” J. Comput. Appl. Math., vol. 122, nos. 1–2, pp. 167–201, 2000. [8] T. Sauer, Y. Xu, and Y. Xuy, “On multivariate Hermite interpolation,” Adv. Comput. Math., vol. 4, no. 1, pp. 1147–1170, 1995. [9] H. Gouraud, “Continuous shading of curved surfaces,” IEEE Trans. Comput., vol. C-20, no. 6, pp. 623–629, Jun. 1971. [10] R. M. Bastos, A. Augusto de Sousa, and F. N. Ferreira, “Reconstruction of illumination functions using bicubic Hermite interpolation,” in Proc. 4th Eurograph. Workshop Rendering, Paris, France, Jul. 1993, pp. 317–326. [11] A. C. Ahlin, “A bivariate generalization of Hermite’s interpolation formula,” Math. Comput., vol. 18, no. 86 pp. 264–273, Apr. 1964. [12] M. M. Chawla and N. Jayarajan, “A generalization of Hermite’s interpolation formula in two variables,” J. Austral. Math. Soc., vol. 18, no. 4, pp. 402–410, Dec. 1974. [13] S. K. Lele, “Compact finite difference schemes with spectral-like resolution,” J. Comput. Phys., vol. 103, no. 1, pp. 16–42, 1992. [14] H. B. Keller and V. Pereyra, “Symbolic generation of finite difference formulas,” Math. Comput., vol. 32, no. 144, pp. 955–971, 1978. [15] W. Burger and M. J. Burge, Digital Image Processing, An Algorithmic Introduction Using Java. New York, NY, USA: Springer-Verlag, 2008. [16] G. Wolberg, Digital Image Warping. Los Alamitos, CA, USA: IEEE Comput. Soc. Press, 1990. [17] A. Belyaev, “On implicit image derivatives and their applications,” in Proc. BMVC, Sep. 2011, pp. 72.1–72.12. [18] A. Belyaev, “Implicit image differentiation and filtering with applications to image sharpening,” SIAM J. Imag. Sci., vol. 6, no. 1, pp. 660–679, 2013. [19] K. K. Delibasis, A. Kechriniotis, and I. Maglogiannis, “On centered and compact signal and image derivatives for feature extraction,” in Artificial Intelligence Applications and Innovations (IFIP Advances in Information and Communication Technology), H. Papadopoulos, A. S. Andreou, L. Lliadis, and L. Maglogiannis, Eds. Berlin, Germany: Springer-Verlag, 2013, pp. 318–327. [20] M. Unser, A. Aldroubi, and M. Eden, “B-spline signal processing. II. Efficiency design and applications,” IEEE Trans. Signal Process., vol. 41, no. 2, pp. 834–848, Feb. 1993. [21] T. Blu, P. Thévenaz, and M. Unser, “Image interpolation and resampling,” in Handbook of Medical Imaging, Processing and Analysis, I. Blankman, Ed. New York, NY, USA: Academic, 2000, pp. 393–423. [22] P. Thévenaz and M. Unser, “Separable least squares decomposition of affine transformations,” in Proc. IEEE ICIP, Santa Barbara, CA, USA, Oct. 1997, no. 200, pp. 26–29. [23] A. Gotchev, K. Egiazarian, and T. Saramaki, “Image interpolation by optimized spline-based kernels,” in Advances in Signal Processing and Applications, J. Astola and L. Yaroslavsky, Eds. Cairo, Egypt: Hindawi Publishing Corp., 2007, pp. 285–335. [24] C. B. Owen and F. Makedon, “Bottleneck-free separable affine image warping,” in Proc. ICIP, Santa Barbara, CA, USA, vol. 1. Oct. 1997, pp. 683–686.

IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 23, NO. 7, JULY 2014

[25] P. Thevenaz, T. Blu, and M. Unser, “Interpolation revisited [medical images application],” IEEE Trans. Med. Imag., vol. 19, no. 7, pp. 739–758, Jul. 2000. [26] K. Delibasis, A. Kechriniotis, N. Assimakis, S. Tassani, and G. Matsopoulos, “Hermite kernels for slice interpolation in medical images,” in Proc. 34th Annu. Int. Conf. Eng. Med. Biol. Soc., 2012, pp. 4369–4373. [27] J. Suckling et al., “The mammographic image analysis Society digital mammogram database,” in Proc. 2nd Int. Workshop Digital Mammography, 1994, pp. 375–378. [28] (2011, Nov. 27). CASIA-IrisV1 [Online]. Available: http://biometrics. idealtest.org [29] M. Unser, P. Thevenaz, and L. Yaroslavsky, “Convolution-based interpolation for fast, high-quality rotation of images,” IEEE Trans. Image Process., vol. 4, no. 10, pp. 1371–1381, Oct. 1995. [30] G. Grevera and J. Udupa, “An objective comparison of 3-D image interpolation methods,” IEEE Trans. Med. Imag., vol. 17, no. 4, pp. 642–652, Aug. 1998. [31] G. Penney, J. Schnabel, D. Rueckert, M. Viergerver, and W. Niessen, “Registration based interpolation,” IEEE Trans. Med. Imag., vol. 23, no. 7, pp. 922–926, Jul. 2004. [32] M. Seppa, “High-quality two-stage resampling for 3-D volumes in medical imaging,” Med. Image Anal., vol. 11, no. 4, pp. 346–360, 2007. [33] J. D. Faires and R. L. Burden, Numerical Methods, Boston, MA, USA: PWS, 1993. [34] R. G. Keys, “Cubic convolution interpolation for digital image processing,” IEEE Trans. Acoust., Speech, Signal Process., vol. ASSP-29, no. 6, pp. 1153–1160, Dec. 1981. [35] D. P. Mitchell and A. N. Netravali, “Reconstruction filters computergraphics,” in Proc. 15 Annu. Conf. Comput. Graph. Interact. Techn., 1988, pp. 221–228. [36] M. Unser, A. Aldroubi, and M. Eden, “B-spline signal processing. I. Theory,” IEEE Trans. Signal Process., vol. 41, no. 2, pp. 821–832, Feb. 1993. [37] T. Blu, P. Thévenaz, and M. Unser, “MOMS: Maximal-order interpolation of minimal support,” IEEE Trans. Signal Process., vol. 10, no. 7, pp. 1069–1080, Jul. 2001.

Konstantinos K. Delibasis received the degree from the Department of Physics, University of Athens, Athens, Greece, in 1990, and the M.Sc. and Ph.D. degrees from the Department of Biomedical Physics and Bioengineering, University of Aberdeen, Aberdeen, U.K., in 1991 and 1995, respectively. He is currently an Assistant Professor with the Department of Computer Science and Biomedical Informatics, University of Thessaly, Volos, Greece. His research interests include processing and analysis of biomedical signal, image, and video.

Aristides Kechriniotis received the degree from the Department of Mathematics, University of Patras, Patras, Greece, in 1975, and the Ph.D. degree from the Department of Engineering Sciences, Polytechnic School, University of Patras, in 2002. He is currently a Visiting Lecturer with the Technological Educational Institute of Sterea Ellada, Lamia, Greece. His research interests include mathematical analysis, applied linear algebra, mathematical physics, and applications to signal processing.