{'fraction_non_alphanumeric': 0.024305555555555556, 'fraction_numerical': 0.00880636070853462, 'mean_word_length': 5.127967930928153, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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.


Introduction
I n 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 Original (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][5][6][7][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][10][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][13][14][15][16][17][18][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][21][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 threedimensional 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 im-age 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]: That assuming small movement and using Taylor series can be developed: Where, u, v, w are components of velocity. I x , I y , I z are spatial and I t 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 spatiotemporal derivations in Equation 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: Where u -k , v -k and w -k denote horizontal, vertical and diagonal neighborhood averages of u k , v k and w k . 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 pro-cessing. 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: 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]:

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 com-mon 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 t 0 and t 1 , we can interpolate images between them using pixel displacements that are obtained by optical flow technique. So, While using Taylor series can be developed as follows: 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 i ) and I(x, t i+1 ): Where t N is the distance between two existed frames N K K is the distance between the first and desired frames N δ  − =   V -1 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 DIRlab 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).
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: 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. 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.
In order to visualize the reconstructed image as a result of four proposed methods, Figure 2 shows the interpolated image at a given slice   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.
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 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. 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).
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.

Discussions
Image reconstruction between two frames of CT image sequences consists of two steps:     ages 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.