A FAST MORPHING-BASED INTERPOLATION FOR MEDICAL IMAGES : APPLICATION TO CONFORMAL RADIOTHERAPY

A method is presented for fast interpolation between medical images. The method is intended for both slice and projective interpolation. It allows offline interpolation between neighboring slices in tomographic data. Spatial correspondence between adjacent images is established using a block matching algorithm. Interpolation of image intensities is then carried out by morphing between the images. The morphing-based method is compared to standard linear interpolation, block-matching-based interpolation and registrationbased interpolation in 3D tomographic data sets. Results show that the proposed method scored similar performance in comparison to registration-based interpolation, and significantly outperforms both linear and block-matching-based interpolation. This method is applied in the context of conformal radiotherapy for online projective interpolation between Digitally Reconstructed Radiographs (DRRs).


INTRODUCTION
Conformal radiotherapy is a cancer treatment protocol that allows to deliver escalating X-ray doses with minimal exposure to surrounding healthy tissues.In order to achieve such goal, the protocol combines state of art technologies allowing the precise anatomical definition and targeting of tumors, as well as the elaboration of dosage goals with input from Computerized Tomography (CT) facilities.
The medical process begins by acquiring a reference CT of a patient, and by computing an irradiation procedure.During several weeks, the patient comes every day for an irradiation session.In order to have a safe and efficient irradiation, the patient must be perfectly positioned, so that the position of the patient at initial CT scan acquisition is reproducible throughout the treatment session.At present, physicians have to accurately estimate patient positioning errors by visual inspection only.However, this is likely to be responsible for limited accuracy and unnecessary loss of time.Therefore, an automated system to characterize and to estimate the displacement of the patient as compared to the reference position has been proposed (Clippe et al., 2003;Sarrut and Clippe, 2001).The method involves the use of portal images generated by Electronic Portal Imaging Devices (EPID) (Fig. 1).Portal Images (PI), which are projective images acquired before the beginning of each treatment session, are compared to the reference CT image through an intensitybased registration between portal images and Digitally Reconstructed Radiographs (DRRs) which are 2D images generated by projection of the 3D CT scan image.Experiments of the automated system leads to a far better accuracy than the manual one, and requires low computation times compatible with a clinical use of the system.However, there are still some steps within the system that need to be further optimized in order to: assure smooth and fast interpolation between the consecutive slices of the 3D CT volume representation in order to provide an isotropic volume image from which the 3D object of interest is reconstructed and thus the irradiation ballistics is accurately calculated.
ensure accurate patient positioning so that the position of the patient at initial CT scan acquisition is reproducible throughout the treatment session by applying the online intensity-based 2D/3D registration between Portal Images and Digitally Reconstructed Radiographs.However, the computing of DRRs can be a time-consuming procedure, especially in the case of intensive use as required for 3D (scanner X) / 2D (PI) image registration.The objective is to optimize the process of DRR generation to fit the time constraints imposed by the clinical application.

Fig. 1. Generation of portal images using Electronic Portal Imaging Devices (EPID)
Image interpolation is of great importance in biomedical visualization and analysis.Threedimensional medical imaging devices usually deliver cross-sectional images of the human body.Piecing together sequences of adjacent images results in a 3D volume representation.The distance between consecutive slices is usually larger than the distance between neighboring pixels within a given slice, the pixels between the slices must therefore be interpolated.
Image interpolation creates a number of new slices between known slices and provides an isotropic volume image from which the 3D object of interest is reconstructed.One example of multidimensional image interpolation is the threedimensional reconstruction of bone geometry and tissue density features.
Interpolation algorithms can broadly be divided into two categories: scene-based (image-based) and object-based (shape-based) algorithms.In scene-based methods, the interpolation is determined only from the image intensities.In object-based methods, some object information extracted from the slices is used for guiding the interpolation process.
The most commonly used scene-based methods are nearest-neighbor and linear interpolations (Goshtasby et al., 1992;Higgins and Ledell, 1994).Although scene-based routines have received a lot of attention (Lehmann et al., 1999;Thévenaz et al., 2000), there has been repeated evidence in the literature of the superior performance of object-based over scene-based interpolation techniques (Barret and Stringham, 1993;Grevera and Udupa, 1998;Herman et al., 1992), especially when the in-plane position of anatomical features is considerably modified between slices.
Object-based interpolation methods have been extended by allowing registration between slices.Probably one of the main registration-based interpolation methods was developed by Penney et al. (2004), where the primary registration step allows spatial correspondence between adjacent slices, then interpolation is carried out between corresponding positions in each slice.
However, despite the efficacy of such technique (Penney et al., 2004), its rather high computational cost can be considered as a serious drawback for clinical applications like conformal radiotherapy where offline-interpolation is not always needed and even if it is, sometimes, it is not possible due to the large number of patients treated at the radiotherapy center.Therefore, our main motivation is to develop an interpolation routine to work in the context of conformal radiotherapy for both slice and projective interpolation.In others words, it allows offline slice interpolation for CT data, and online projective interpolation between Digitally Reconstructed Radiographs (DRRs).The method has to be flexible, simple to implement and to evolve, accurate but without the disadvantage of a more expensive computation cost.
In a previous investigation, we presented some early experiments of the method that showed its usefulness for projective interpolation between DRRs (Atoui et al., 2004).In this paper, we present and assessed the method for slice interpolation of 3D CT data.The method shares some of its philosophy from registration-based methods but differs from them by using a less costly mapping routine: correspondence is established using a block matching process.It also uses an image morphing algorithm to interpolate between adjacent slices.This morphing-based interpolation will be compared to standard linear interpolation, blockmatching-based interpolation and registration-based interpolation in tomographic data sets.

Interpolation between consecutive slices involves a two-step formalism:
-The first step is matching, which is based on establishing a correspondence between source and target images via block matching.
-The second step is interpolation, which uses the correspondence obtained in the first step, and generates the in-between image through image morphing.

MATCHING STEP
In medical images, the pixels corresponding to soft tissue contours and bone structures have high gradient values and thus can be used as anatomical landmarks.Using a block matching process, we can establish the correspondence between two adjacent images, based on the differences of pixel intensity and gradient.
The block matching algorithm operates as follows: -Each image frame is divided into a fixed number blocks.
-For each block in the source image, a search is made over a limited area in the target image to obtain the best matching block with the least prediction error.
-A block is considered found if it matches another block in the second image with a correlation higher than a threshold T .
-The block matching algorithm results in two meshes of control points: coordinates of the upper right-hand corners of the blocks.The source image is associated with a mesh M S which specifies the coordinates of control points.A second mesh M D specifies their corresponding positions in the target image.
The efficacy of the matching depends on four parameters: the block size B, the search region R, the threshold T and the similarity criterion S. Values of B = 16, R = 8 and T = 0.85 were found to provide a reasonable compromise for the sequences tested.Similarity was defined from the intensities and gradients of the blocks.In other words, matching a block in the target image to a block in the source image is a normalized function of the intensity, the gradient magnitude, and the gradient direction of the block as follows: where I(x,y), M(x,y) and D(x,y) are, respectively, the intensity, the gradient magnitude, and the gradient direction of the block (x, y) in the source image (note that x and y are the coordinates of the upper righthand corner of the block, and the similarity measure is applied to the whole block); and I (x , y ) M (x , y ) D (x , y ) are, respectively, the intensity, the gradient magnitude and the gradient direction of the block (x , y ) in the target image.

INTERPOLATION STEP
Morphing between two images involves the deformation of reference images towards an intermediate position.It consists in generating two intermediate images, followed by the application of a color blending technique (cross dissolve, which is equivalent to linear interpolation) to compute the final artificial image.
Image morphing provides flexibility for accommodating local transformations between images and makes it possible to generate continuous deformations from source image to destination image.
Numerous morphing techniques (mesh warping, field morphing, energy minimization, free-form deformations, etc) (Wolberg, 1998) are available for the generation of artificial intermediate images.Generally, all these techniques follow similar procedures, including establishing correspondences between images (via points, segments, curves, etc), calculating mapping functions, controlling transition and finally pixel interpolation (cross-dissolve).However, a simple, rapid, efficient and reliable morphing technique was desirable to investigate the feasibility of our approach.For this initial platform, we have chosen to use the 2-pass mesh warping algorithm (Wolberg, 1990).
The mesh-warping algorithm relates features with the mesh of control points provided by the matching process as follows: -The algorithm accepts a source image (Fig. 2d), a destination image (Fig. 2f) and two arrays of coordinates.The first array, M S , specifies the coordinates of control points in the source image (Fig. 2a).The second array, M D , specifies their corresponding positions in the destination image (Fig. 2i).Both M S and M D must have the same dimensions in order to establish a one-to-one correspondence.
-Then two images are processed through a 2-pass warping with 2 output intermediate images I 1 (Fig. 2c) and I 2 (Fig. 2g).
• The first pass is responsible for resampling each row independently.It maps all initial image points coordinates (u, v) to their (x, v) coordinates in the intermediate image, thereby positioning each input point into its proper output column.
• The second pass then resamples each column in the intermediate image, mapping every (x, v) point to its final (x, y) position in I 1 and I 2 .Cubic splines (Catmull-Rom) are used to interpolate new pixel locations between warped mesh vertices.Catmull-Rom spline interpolation allows for smooth interpolation between two different meshes to determine the correspondence between pixels.Fant's algorithm (Fant, 1986) is used to resample the pixel color according to the spatial mapping provided by Catmull-Rom spline interpolation.
-After the landmarks of the source and target images have been aligned through (I 1 and I 2 ), they are cross-dissolved to generate a final in-between image I (Fig. 2e).In other words, the intensity value of the pixel u(x, y) in the middle image I is the result of the linear interpolation between u 1 (x, y) of I 1 and u 2 (x, y) of I 2 .

RESULTS
The purpose of this section is to provide a twofold analysis of the performance of the proposed method: firstly, with a quantitative error analysis and comparison with linear and block-matching-based interpolation; secondly, with a visual evaluation of the quality of interpolated images.
To determine how well the block matching is doing, we conducted some experiments in terms of landmark localization.The scenario was to consider two consecutive slices in a 3D CT image, place several anatomical landmark points on the first slice and the corresponding points (determined visually) on the second slice.Then the block matching algorithm used the mesh of points of the first slice and determined the corresponding mesh in the second slice.Finally, we calculated the distance between the reference points and the points obtained with the block matching algorithm.As these experiments were not complete, A series of tests were then conducted to evaluate the quality of the images generated by morphing as compared to reference images.For this study, 9 3D CT abdomen and 2 3D CT lungs images were available.Image dimensions and voxel sizes are summarized in Table 1.All images were acquired from scanners performed at Léon Bérard anti-cancer center.
We reused the framework for error analysis proposed by Grevera and Udupa (1996).Our approach was tested against three interpolation techniques: linear interpolation, block-matching-based interpolation and registration-based interpolation.
As with our method, the block-matching-based interpolation has a matching step using block matching.However, in the interpolation step, the morphing is replaced by simple linear interpolation between the corresponding blocks of pixels in source and destination slices.
As for the registration-based interpolation, we used a home-made implementation of a 2D deformable registration algorithm (based on Sarrut et al., 2006) to obtain a dense deformation field.Our method is slightly different of Penney et al. (2004) registrationbased interpolation because 1) we used Sum of Squared Differences (SSD) instead of Normalized Mutual Information as similarity criterion, and 2) we used Gaussian regularization instead of B-splines parametrization.Moreover, the resulting vectors field is linearly interpolated from a 3×3 mm 2 grid although, in Penney et al. (2004), it was interpolated with Cubic B-splines from 20 × 20 mm 2 and 10 × 10 mm 2 grid resolution.Then, slices were interpolated using the bilinear interpolation proposed in Penney et al. (2004).
The assessment implies using sequences of three images and interpolating an in-between image between the first and the third one.The interpolated slice is then compared with the second image, using three error measures: -Mean-squared difference (MSD): where N S is the number of slices being interpolated and N P the number of pixels per interpolated slice.I i (u) and I i (u) m are respectively the original density value and the density estimation obtained using the interpolation method.
-Number of sites of disagreement (NSD): where -Largest difference (LD): In addition to these measurements, we used another measure from Grevera and Udupa (1996), called statistical relevance r.This parameter shows the percentage difference between two interpolation techniques (i.e.m 1 and m 2 ) for a given error measure (i.e.err) as depicted in the following equation: Positive values of r m 1 ,m 2 indicate that the first interpolation technique performs better than the second interpolation technique.Negative values show the opposite effect.
One final test to estimate the matching quality for the interpolated slice in comparison to the original slice, for linear, block-matching and morphing-based interpolation methods, is the mean sum of absolute difference (MSAD) as follows: Error measurements provided by Eqs.1-4 and 6 for linear interpolation (LIN), block-matching-based interpolation (BLIN), registration-based interpolation (REG) and morphing-based interpolation (MORPH) are summarized in Table 2 and Table 3.The values are computed from the 3D CT images as follows: The results in Table 2 and Table 3 show that, in average, our approach gives better estimates of original density values than both linear and block-matchingbased interpolation.Moreover, similar performances in terms of the used metrics were achieved between registration-based interpolation and the proposed morphed-based method.
We are aware that only limited datasets were used to perform the experiments (in particular only CT images were used) and that larger experiments should be performed to compare registration-based and morphed-based interpolation (for example with MR images).However, with more than 800 interpolated slices, computed criteria were very similar between the two methods and significantly different from the linear (LIN) or BLIN ones Fig. 3 illustrates the differences between the interpolation techniques.The problem with linear interpolation is that if features do not line up exactly, interpolation results in a blurred image.In morphingbased interpolation, features from source and target images have already been aligned using block matching.Also, comparing to block-matching-based interpolation, the morphing step in our method allows for smoother transition between images.Visually, morphing-based interpolated images are therefore much similar to the original images.

COMPUTATIONAL COST
When comparing the computation times between our morphing-based interpolation approach, and the registration-based method, we noticed that interpolation between two 512 × 512 slices using our approach took about 15 sec in comparison to about 1 minute for the registration-based interpolation, as measured on a Pentium IV, 2 GHz computer.
However, a quality-based comparison using the same CT database is needed to make sure that our method delivers as accurate images as the registrationbased methods.
The computation time of our method is clearly dominated by the matching step.In other words, the total timing strongly depends on the block matching parameters (block size, search region . . .).One of the possible speedup options for reducing the computation time would be the use of an optimized block matching algorithm as in Salari and Li (1995).

EXISTING ROUTINES
When compared to both shape-based and linear interpolation, registration-based interpolation allows for a statistically better synthesized images Penney et al. (2004).However, while in some cases, images can be interpolated off-line, in other cases, online interpolation is a necessity.In conformal radiotherapy context, to determine the patient positioning before radiation, the whole process: from acquiring the portal images, to conducting the registration between PI and DRRs, to adjusting the patient positioning has to be done in few minutes.Therefore, a simple, easy to deploy, flexible and fast reconstruction routine is preferred especially in a clinical context as conformal radiotherapy.

FUTURE WORKS
In the context of conformal radiotherapy, the success of an automatic repositioning system is closely related not only the accuracy of repositioning the patient, but also to computation time needed to conduct the process so it can compatible with a clinical use of the system.
Basically, to calculate the positioning of the patient at the beginning of the treatment session in order to adapt the treatment process, a solution is to apply an online intensity-based 2D/3D registration (Clippe et al., 2003, Fig. 4) between Portal Images representative of the patient positioning at the treatment session, and Digitally Reconstructed DRRs are computed from the patient CT data with volume rendering (Li and Miguet, 1992;Lacroute, 1995), thus permitting to produce a virtual patient model.DRRs being digitalized, they can easily be managed by means of computer systems.
However, the computing of images can be a time-consuming procedure, especially in the case of intensive use as required for 3D (scanner X) / 2D (PI) image registration.For example, to cover a displacement interval of [−5 • , +5 • ] and a half-degree precision, it is necessary to generate 900 DRRs per view per patient.It is therefore crucial to devise a solution for the optimization of the generation process.
One possible solution would be to generate a set of DRRs off-line representative of a large displacement interval but with an average precision, then at the beginning of every treatment to conduct the following procedure: conduct a first 3D/2D registration between the DRRs and the PI images.The registration will indicate the area of interest, in terms of best match between the portal image and the set of DRRs.
locally, at the area of interest, generate artificial inbetween DRRs via our interpolation method.This provides a high degree of accuracy locally with a reasonable number of DRRs to be generated offline and online comparing to an exponentially growing amount of costly DRRs that had to be computed offline to allow a high degree of precision for the whole displacement interval to apply a one step registration.
Therefore, apart the usefulness of our method in slice interpolation, a second use is projective interpolation where image morphing allows the creation of new artificial in-between images from relatively less similar images like DRRs.In fact, in a previous investigation, we presented some early experiments of the method that showed its usefulness for projective interpolation between DRRs (Atoui et al., 2004).

CONCLUSION
This paper introduces a simple, fast, and flexible morphing-based interpolation algorithm of grey-scale tomographic data by 2D deformation.We compared the performance of our method to linear, blockmatching and registration based interpolation methods for 3D CT images and concluded that our method scored similar performance to registration-based interpolation, and significantly outperforms both linear and block-matching-based interpolation.
The main strength of our method lies in the simplicity of implementation, using well established, easily understandable and optimized techniques (block matching, mesh warping, . . .).Another advantage of the method is its usefulness in both slice and projective interpolation.
Studies are ongoing to compare the method with state-of-the-art routines (registration-based interpolation. . .), to improve our method by using a more accurate and more optimized matching process, and to assess the method within the overall scenario of patient positioning in conformal radiotherapy.

Fig. 2 .
Fig. 2. Illustration of the 2-pass mesh warping.The algorithm allows the deformation of two reference images towards an intermediate selection, thus generating two intermediate images, followed by the application of a cross dissolve to compute the final in-between artificial image.

Fig. 4 .
Fig. 4. The patient displacement error is calculated by comparing 2 PIs with the 3D CT scan.

Table 1 .
Summary of the image dimensions and voxel sizes for the two data set types

Table 3 .
Error measurements for the linear interpolation (LIN), the block-matching-based interpolation (BLIN), the registration-based interpolation (REG), and the morphing based interpolation (MORPH).A dash indicates that the difference between the two techniques was not statistically significant (t-test ≤ 0.01)conduct a second registration between the portal image and the DRRs within the area of interest after adding the newly generated in-between DRRs through morphing.