A Study on Robustness of Various Deformable Image Registration Algorithms on Image Reconstruction Using 4DCT Thoracic Images.

BACKGROUND: Medical image interpolation is recently introduced as a helpful tool to obtain further information via initial available images taken by tomography systems. To do this, deformable image registration algorithms are mainly utilized to perform image interpolation using tomography images.
MATERIALS AND METHODS: In this work, 4DCT thoracic images of five real patients provided by DIR-lab group were utilized. Four implemented registration algorithms as 1) Original Horn-Schunck, 2) Inverse consistent Horn-Schunck, 3) Original Demons and 4) Fast Demons were implemented by means of DIRART software packages. Then, the calculated vector fields are processed to reconstruct 4DCT images at any desired time using optical flow based on interpolation method. As a comparative study, the accuracy of interpolated image obtained by each strategy is measured by calculating mean square error between the interpolated image and real middle image as ground truth dataset.
RESULTS: Final results represent the ability to accomplish image interpolation among given two-paired images. Among them, Inverse Consistent Horn-Schunck algorithm has the best performance to reconstruct interpolated image with the highest accuracy while Demons method had the worst performance.
CONCLUSION: Since image interpolation is affected by increasing the distance between two given available images, the performance accuracy of four different registration algorithms is investigated concerning this issue. As a result, Inverse Consistent Horn-Schunck does not essentially have the best performance especially in facing large displacements happened due to distance increment. Copyright: © Shiraz University of Medical Sciences.

Introduction

In recent years among medical image processing algorithms, image interpolation has been taken into account as a powerful tool due to several applications in both diagnostic and therapeutic fields. When taken CT images from a real patient are not the same quality as desired for clinical diagnosis, image interpolation may be helpful to generate further information for possible lesions. Furthermore, this technique is therapeutically very useful in Image Guided Radiotherapy (IGRT) where additional image information is necessary for precise target localization in order to enhance treatment quality [1,2]. At IGRT, targeting accuracy will be a crucial issue when tumors located in thorax region of body move mainly due to breathing phenomenon [3]. Moreover, by using image interpolation, the patient is kept away from additional dose of re-scanning for getting new images considering ALARA principle.

Several strategies have been proposed in different studies on generation a set of spatial interpolated images at any arbitrary time ranging from shape-based to intensity-based interpolation algorithms [4-8]. Additionally, image registration tool can also be considered as an alternative approach to interpolate between two images.

While taking CT images from dynamic organs located in thorax, motion information of all including objects is missing between two consequent images [9-11]. In order to achieve such information, image registration technique has been proposed in different features. Image registration gives a unique pattern including displacement field between two existing images. A range of image registration techniques, including rigid and non-rigid registration, has been developed in order to find the information of mismatched or deformed organs between such images [12-19]. Conceptually, the continuous displacement field obtained by image registration technique may be optimal for reconstructing interpolated images at any time. In this study, the robustness of commercially available image registration techniques will be taken into account for image reconstruction due to intrinsic properties of each technique where no comparative study has been performed before. For this purpose, we focused on two highly performing intensity-based deformable registration algorithms for lung CT images known as Demons (Thirion 1998) and the original H.S. optical flow (Horn & Schunck 1981) algorithms [20-22]. Two Horn-Schunck (Original HS, Inverse Consistent HS) and two Demons (Original Demons, Fast Demons) image registrations were evaluated by means of DIRART algorithm. DIRART software package is a free software dedicated for deformable image registration and adaptive radiotherapy research developed by Dr. Deshan Yang [23].

Since image interpolation is highly affected by the distance between two given images, the performance accuracy of different registration algorithms was investigated by increasing the distance between two source and target images. In other words, we are interested in assessing “how well deformable registration algorithm can perform interpolation by increasing the spatial distance between two given images”. To address this issue, several distances were tested and compared with each other.

Thoracic images of five real patients provided by DIR-lab [24] were used in this work. Dataset of each patient consisted of three-dimensional CT data at 10 breathing phases. DIR-lab images including a large range of reference samples with different spatial distributions have been proposed for investigating the accuracy of DIR performance. In order to evaluate the performance of each registration strategy at image reconstruction, the real image located in the middle part of two source and target images was chosen as ground truth data. After implementation of our code, image output is compared with real middle image and the differences will be discussed using conventional mathematical approach.

The final analyzed results represent that Inverse Consistent Horn-Schunck algorithm with the least error has the best performance in reconstructing interpolated images at lower computational time making it very promising in clinical practice.

Material and Methods

This work represents a comparative study that analyzes the performance of deformable medical image registration algorithms using 4DCT, quantitatively. In general, medical image registration algorithms present a map that illustrates the correspondence of different features between two medical images taken by medical imaging systems. When a change on spatial shape of organ or volume of interest happens in thorax region, a complex non-rigid algorithm known as deformable image registration (DIR) can be sufficient for an accurate aligning correction against simple rigid approximations. DIR technique is divided into two intensity-based and featured-based image registration techniques. Intensity-based method that works by optical-flow-like methods is a fully automatic algorithm using the intensity distribution of the two images for statistical measures of similarities [25,26]. In this section, the utilized intensity-based deformable registrations consist of Horn-Schunk, and Demons algorithms will be briefly considered to provide data required for image reconstruction performing.

Optical Flow

The goal of optical flow methods is to solve motion equation acquired between two frames which are taken at times t and t+Δt at every voxel position. For a voxel with intensity I(x,y,z,t) that has moved by Δx, Δy, Δz and Δt between two image frames, motion equation can be given [27]:

I(x,y,z,t)=I(x+Δx,y+Δy, z+Δz, t+Δt) (1)

That assuming small movement and using Taylor series can be developed:

I (2)

Where, u, v, w are components of velocity. I are spatial and It is temporal image derivative. This is known as the aperture problem of the optical flow algorithms that cannot be solved as such.

Many methods have been suggested to solve this equation [28]. Differential methods based on partial derivatives of image signal are more applicable.

Horn-Schunck Method

The Horn-Schunk (HS) algorithm proposed by Horn and Schunck [29] is a global optical flow method with an accurate performance. HS algorithm is based on differential technique that uses gradient constraint with a global smoothness to obtain velocity field. This algorithm consists of two steps.

At first, the gradient constraint with a global smoothness is used in order to estimate spatio-temporal derivations in Equation 3.

(3)

The first term of this equation shows error in brightness constancy and the second term represents global smoothness, so the whole equation illustrates errors or distortions in flow. HS algorithm tries to minimize distortions in flow and prefers solutions which show more smoothness; therefore, in the next step, sum of the errors are minimized by solving iterative equation (3) in order to obtain final motion vector. Finally, components of velocity are given by:

(4)

(5)

(6)

Where u and w denote horizontal, vertical and diagonal neighborhood averages of u and w. The superscript k+1 denotes the next iteration, which is to be calculated and k is the iteration number. The smoothness weight (α) plays an important role where the brightness gradient is small, for which the suitable value should be determined.

Demons Method

The concept of demons was introduced by Maxwell to illustrate a paradox of thermodynamics. And then, Maxwell’s demon was adapted by Thirion [30] to use in image processing. In this method, diffusion algorithm was utilized to align two moving and reference images. Assuming I and J as reference and moving images respectively, the aim is J deformation to be similar to I as much as possible. This technique uses gradients of moving and reference images which determine the direction of each voxel. The deformation field is smoothed by a Gaussian filter, and iteratively is used to transform the moving image, and register on to the reference image. Finally, the displacement field consists of individual vectors corresponding to each voxel. The moving image is iteratively deformed by applying a displacement vector dr = (dx, dy, dz) to each voxel as:

(7)

Where J(n) and I(n) are the intensity of the moving and reference image at the n-th iteration; J(0) and I(0) are the original intensity of the moving and the reference image.

Fast Demons Method

Equation (7) only uses gradient information from a reference image to determine the demon force, and it can cause problems when the gradient of the reference image is small [31]. This problem may be corrected using the gradient of the iteratively updated moving image [32]:

(8)

Inverse Consistent Deformable Image Registration

Deformable image registration is called inverse consistent if there is no difference between given source and target images. Considering U and V as forward and backward transformations between I and J respectively, following equations are applied:

I◦U=J and J◦V=I (9)

Inverse-Consistency is presented as common and more accurate registration algorithm. Dr. Deshan Yang proposed a new algorithm that is more accurate and faster than previous suggested inverse consistent algorithms [33]. Both images register towards each other until both deformed images are matched and register correctly. In each pass, images are deformed with the delta motion field that is acquired using minimizing a symmetric optical flow cost function on positive and negative directions.

In this work, Inverse consistent HS, Original HS, Fast demons and Original demons are investigated by means of DIRART software package. The value of smoothness during iteration, α2 for HS algorithm and the Gaussian low-pass filter window size for Demons were set to 3 in real voxel sizes. Also, we used the max filter to compute the image intensity during image down sampling in order to obtain better results [23].

Interpolation Method

Interpolation is the process of estimating new values within the range of known values, being commonly used in medical image processing. Spatial interpolation is the process of estimating the value of unknown points within the object’s area using existing points while temporal interpolation is the estimation of the value of an object at a time point using data from nearby time points. The proposed method by Jan Ehrhardt [4] is a temporal interpolation method used in this study.

Optical Flow-based Interpolation Method

Given two images at time t0 and t1, we can interpolate images between them using pixel displacements that are obtained by optical flow technique. So,

I(x(t), t, (10)

While using Taylor series can be developed as follows:

I(x(t),t (11)

But in general, the intensity conservation assumption might not be fulfilled, and structures may appear or disappear between two time steps. In this work, we used a weighted average between the corresponding voxel intensities in the adjacent time frames I(x, t and I(x, t:

(12)

Where

V cannot be computed directly. It can be obtained using methods such as gradient descent, Gauss-Newton or other more stable iterative methods. In this study, Chen [34] method for four deformable image registration technique was taken into account.

4DCT Patient Database Properties

Dataset used in this work includes thoracic 4DCT images of five patients taken at the University of Texas M. D. Anderson Cancer Center in Houston. This dataset assessed by DIR-lab group includes 128 slices with 2.5mm slice thickness acquiring with a General Electric CT scanner (GE Medical Systems, Waukesha, WI). Each slice is a 2D image with different dimensions and voxel sizes for each patient uniquely (Table 1).

Table 1

CT Images Characteristics

Patient	Image Dimension	Voxel Size(mm)
1	256 x 256 x 94	0.97 x 0.97 x 2.5
2	256 x 256 x 112	1.16 x 1.16 x 2.5
3	256 x 256 x 104	1.15 x 1.15 x 2.5
4	256 x 256 x 99	1.13 x 1.13 x 2.5
5	256 x 256 x 106	1.10 x 1.10 x 2.5

Total 4DCT data of each patient consists of a set of 3D CT images of ten points during whole breathing cycle. For example, CT1 represents a set of three dimensional images taken at maximum inhale phase of breathing, while CT6 corresponds to the maximum exhale phase. Therefore, CT1 to CT5 correspond to the inhale phase and CT6 to CT10 correspond to the exhale phase.

Results

In order to assess the robustness of available deformable image registration methods chosen in this work to perform image reconstruction, the generated image has been compared with real middle image as benchmark and the differences are presented by means of Mean Square Error (MSE) quantitatively and also visual difference images. For this aim, CT2 and CT4 were selected as existing images and CT3 was assumed to be used versus images reconstructed via proposed image registration algorithms. MSE for actual (I(x,y,z,t)) and interpolated (J(x,y,z,t)) images with M×N×X size is defined as:

(13) Figure 1 shows the calculated mean square error between a set of 3D constructed images (depending on the number of slices at each database) and three dimensional middle image databases (CT3). For example, 94 images will be reconstructed and compared with 94 real images of CT3 dataset of patient No.1. This value is 648mm, 656mm, 1964mm and 724mm using Original HS, Inverse consistent HS, Original Demons and fast Demons algorithms, respectively at reconstructing chosen 2D paired images among 3D image dataset of patient 1.

Figure1

MSE between Interpolated and Real Middle Images Using four DIR Algorithms over five Patients

Table 2 reports the MSE of reconstructed and real images at middle slice for the first patient. As seen in this table, Inverse consistent HS algorithm demonstrates the best performance in reconstructing the image.

Table 2

Mean Square Errors of Reconstructed and Real Images at Middle Slice for the First Patient

DIR Methods	MSE
Original HS	648.41
Inverse Consistent HS	656.1
Original Demons	1964.2
Fast Demons	724.4

In order to visualize the reconstructed image as a result of four proposed methods, Figure 2 shows the interpolated image at a given slice belonging to image database by Original HS (a), Inverse consistent HS (b), Original Demons (c) and Fast Demons (d) registration algorithms in comparison with real middle image as a benchmark for the first patient (e).

Figure2

Reconstructed Images Generated by Original HS (a), Inverse Consistent HS (b), Original Demons (c) and Fast Demons (d) Algorithms against Real Middle Image (e)

Figure 3 shows the differences between interpolated images and real middle image at the same slice of a given patient as different image. As seen in this figure (Figure 3), less difference resulting in better matching between two interpolated and real image is given by Inverse Consistent HS algorithm (Figure 3b). In ideal case where there are no differences, a uniform gray image without any contrast must be resulted.

Figure3

Differential Images between Interpolated Images Generated Using Original HS (a), Inverse Consistent HS (b), Original Demons (c) and Fast Demons (d) Registration Algorithms and Real Middle Image (CT3)

To assess the performance accuracy of different registration algorithms by increasing the distance between two given source and target image datasets, 3D CT No. 5 was considered as ground truth image and two 3D CT image databases before and after CT No. 5 were selected to represent distance increment: CT4 and CT6 with less distance, CT3 and CT8 with median distance, CT1 and CT9 with large distance. It should be noted that increasing distance may disturb the smoothness degree of displacement.

Table 3 shows the effect of increasing the distance between two sets of 3D CT images for patient 1. As shown, the Inverse consistent HS algorithm has better performance in small distances. However, this algorithm is not able to work as well by increasing distance. The least MSE belongs to Original HS for low smoothness degree while the distance between two images is increasing.

Table 3

MSE of Reconstructed Image by Four Proposed Methods vs. CT5 Using Two Different CT Data Pairs

	MSE(mm) at different DIR
Two Selected 3D CT Images	Original HS	Inverse Consistent HS	Fast Demons	Original Demons
CT4_CT6	712	710	1561	2053
CT3_CT8	4094	4056	4472	4857
CT1_CT9	8539	8544	8751	8986

Table 4 represents the MSE between two assumed 2D images among 3D CT data of patient 1 considering the same calculations used above for assessing distance increasing effect between two 2D source and target images.

Table 4

MSE among Two 2D Reconstructed and True Image by Four Proposed Methods vs. 2D CT No.5 for Patient 1

	MSE(mm) at different DIR
Two Selected 2D CT Images	Original HS	Inverse Consistent HS	Fast Demons	Original Demons
CT4_CT6	240	236	244	1698
CT3_CT8	3815.5	3815.5	3850	4560
CT1_CT9	6432	6313	8407	9894

Figure 4 shows reconstructed images generated by Inverse consistent HS between CT4-CT6 (a), CT3-CT8 (b) and CT1-CT9 (c) and compared with a given real slice of CT5 (d).

Figure4

Reconstructed Images versus Corresponding Real Image (d) Using Two Image Pairs as CT4-CT6 (a), CT3-CT8 (b), CT1-CT9 (c)

Different images emerged between reconstructed images and real image are shown in Figure 5 to give a better depiction of performance accuracy regarding distance increment. As resulted in this Figure (Figure 5a) and (Table 3), the reconstructed image with less difference is derived while the distance is in minimum value.

Figure5

Different Images between Actual CT5 Image and Interpolated Images Emerged between CT4-CT6 (Figure 5-a), CT3-CT8 (Figure 5-b) and CT1-CT9 (Figure 5-c)

Discussions

Image reconstruction between two frames of CT image sequences consists of two steps: first, finding the displacement field between two given images and then, the optical flow-based interpolation is used to generate an image at the desired time. Therefore, the registration algorithms play an important role in image interpolation with high accuracy. Since the most common algorithms to implement deformable image registration for lung CT images are Horn-Schunck and Demons algorithms, we assessed the abilities of two Horn-Schunck and two Demons algorithms in image reconstruction by means of commercial DIRART software package developed by MATLAB. By comparing the performance of Original Horn-schunck, Inverse consistent Horn-schunck, Original Demons and Fast Demons algorithms, Horn-Schunck based algorithms may be optimal against Demons-based algorithms in the reconstruction of thoracic images.

Both HS and Demons algorithms use gradient to determine the direction of each voxel but their strategy at smoothing deformation field is uniquely different due to their intrinsic properties resulting in reconstructed images with different accuracies in details. Another difference between these two algorithms is their computational time. In each iteration, Demons algorithm run times are longer than HS algorithms due to re-sampling deformed image, while the number of re-sampling steps in HS algorithm is few to make it promising in practical cases.

Between two presented HS algorithms, the better performance belongs to Inverse consistent HS algorithm in generating virtual images due to generating symmetric registration.

It should also be noted that Fast Demons algorithm had better operation in comparison with Original Demons due to using the gradient of the iteratively updated moving image.

Moreover in this study, the effect of increasing distance between two source and target images on the accuracy of image reconstruction was taken into account. Final analyzed results showed that Inverse consistent HS registration algorithm proved to be the best image reconstruction accuracy for small distances.

The performance of commercially available HS and demons algorithms in image reconstruction were investigated and tested in this study using database of real patients. Results showed that Inverse consistent HS approach yields the best estimation of desired image against Original Demons. However, higher levels of accuracy can even be obtained modifying optimizing parameters (e.g. applied optimal filters) considering model simplicity to avoid long computational time.

A future avenue for research can be the parameter optimization of HS algorithms to find a more accurate image reconstruction.

Conclusion

In this work, a comparative study was done to assess the role of four deformable image registration algorithms for generating virtual interpolated image between two consequent tomography images. To do this thoracic images of five real patients at 10 breathing phases were utilized. The virtual interpolated image generated by deformable registration algorithms were compared with real image located in the middle part of two source and target images, as ground truth data. The final obtained results show that Inverse Consistent Horn-Schunck algorithm with the least error has the best performance in reconstructing interpolated images. Moreover, lower computational time of this algorithm makes it very promising for clinical practice.