A fast 4D cone beam CT reconstruction method based on the OSC-TV algorithm

Abstract

BACKGROUND:

Four-dimensional cone beam computed tomography allows for temporally resolved imaging with useful applications in radiotherapy, but raises particular challenges in terms of image quality and computation time.

OBJECTIVE:

The purpose of this work is to develop a fast and accurate 4D algorithm by adapting a GPU-accelerated ordered subsets convex algorithm (OSC), combined with the total variation minimization regularization technique (TV).

METHODS:

Different initialization schemes were studied to adapt the OSC-TV algorithm to 4D reconstruction: each respiratory phase was initialized either with a 3D reconstruction or a blank image. Reconstruction algorithms were tested on a dynamic numerical phantom and on a clinical dataset. 4D iterations were implemented for a cluster of 8 GPUs.

RESULTS:

All developed methods allowed for an adequate visualization of the respiratory movement and compared favorably to the McKinnon-Bates and adaptive steepest descent projection onto convex sets algorithms, while the 4D reconstructions initialized from a prior 3D reconstruction led to better overall image quality.

CONCLUSION:

The most suitable adaptation of OSC-TV to 4D CBCT was found to be a combination of a prior FDK reconstruction and a 4D OSC-TV reconstruction with a reconstruction time of 4.5 minutes. This relatively short reconstruction time could facilitate a clinical use.

Keywords

4D CBCT iterative reconstruction respiratory motion multi-GPU

1 Introduction

External beam radiotherapy in the thoracic region presents the particular problematic of a substantial target displacement during the treatment due to respiratory motion [1]. To reduce the irradiation of healthy tissues, image-guided radiation therapy (IGRT) should consider this temporal displacement by using a 4D representation of the patient’s body. While computed tomography (CT) -based 4D imaging is now routinely used in the clinic to better define the treatment target, 4D imaging is not yet commonly used at the time of treatment for IGRT. Being acquired over a relatively long period, cone beam computed tomography (CBCT) projections contain information about this displacement on a day-to-day basis. A 4D reconstruction algorithm aims to recover this information and use it to improve the general quality of the reconstructed image and, eventually, the accuracy of IGRT.

A straightforward approach to achieve a time correlated CBCT would be to simply sort projections according to the respiratory phase during which they were acquired before reconstructing each phase independently. However, the low number of projections from each phase and their uneven spacing would lead to major artifacts, particularly if the reconstruction is performed by a conventional filtered backprojection algorithm such as FDK [7]. To reduce these artifacts, the McKinnon-Bates (MKB) [17] algorithm offers to correct the 3D FDK image using phase reconstructions of the difference between the initial raw data and the forward projected data.

The vast field of iterative reconstruction, enabling the use of physics modeling and a priori knowledge of image properties, offers many promising possibilities in 4D imaging. Total variation (TV) regularization [27] is particularly suitable to partially correct the streaking artifacts and notably improve the global image quality [21, 22]. Considering that many image regions remain static during respiratory movement, the use of information from all projections to reconstruct each phase can also improve image quality. To this end, different strategies have been explored, such as optical flow based registration [6] or auto-adaptive phase correlation algorithm [3]. The use of a prior 3D reconstruction had also been considered to constrain the convergence of an iterative reconstruction algorithm, using for example compressed sensing [5, 14]. Wang et al. have used a planning CT and deformation vector fields (DVF) to obtain high-quality CBCT [31], an approach that was later extended to enable simultaneous motion estimation and image reconstruction [32], notably with a mesh-based organ motion modeling [36]. The adaptive steepest descent projection onto convex sets (ASD-POCS) algorithm [28] is another popular method in 4D CBCT which uses compressed sensing.

Beyond reconstruction algorithms alone, compound task-driven image acquisition and reconstruction have been recently proposed by Gang et al. [8]. This method optimizes the acquisition parameters of each projection view and its subsequent filtering, while relying on a mathematical model of task-based image quality criteria. This approach yields promising results, while heavily relying on exceptional flexibility of the imaging equipment (C-arm CT with custom scanning protocol), which is not always clinically accessible.

The present work studies a novel method designed to improve image quality and ease of clinical implementation of 4D CBCT by combining different approaches. To be easy to use in clinic, a 4D CBCT algorithm should use a standard set of projections (no slower rotating gantry or larger number of projections acquired), and should not require additional steps from the operators, such as manual delimitation of mobile regions. Reconstruction should be as fast as possible, ideally in real-time, to allow for treatment plan correction while the patient is in the treatment position. Real-time 4D CBCT reconstruction appears unlikely in the short term, and typical reconstruction times are reported in the 5 min – 30 min interval [18 , 35]. The current paper is focused on offering reasonable reconstruction times of a few minutes.

The selected reconstruction algorithm, OSC-TV [16], combines a modified ordered subsets convex algorithm (OSC) [11] and a total variation minimization regularization technique. The OSC algorithm is particularly fast, which would facilitate a clinical use. The TV regularization reduces noise and undersampling artifacts, which could help to obtain sufficient image quality without a slower gantry rotation or a multiple revolutions protocol. To adapt the OSC-TV algorithm to 4D imaging, the introduction of knowledge from all projections through the initialization of each phase’s reconstruction by a prior 3D reconstruction was studied. This alternative initialization scheme is compared to a standard initialization by a blank image using reconstructions of a numerical phantom and of clinical datasets. The different methods of adapting OSC-TV to 4D imaging are compared to a filtered backprojection approach, MKB, and an iterative one, ASD-POCS. In order to reduce computation time, the parallelization of the phases’ reconstruction using several GPUs was also studied. Two main objectives are pursued in this study: the adaptation of the OSC-TV algorithm to 4D CBCT to rapidly obtain a sufficient image quality and the comparison of the aforementioned initialization schemes.

2 Materials and methods

2.1 OSC-TV algorithm

The selected reconstruction method, the OSC-TV technique [16], is an expectation-maximization algorithm based on Poisson photon statistics with promising properties for a fast 4D CBCT reconstruction of conventionally acquired 3D CBCT. Such an iterative algorithm should be able to reconstruct an image of better quality than an analytic approach like MKB. For 3D CBCT imaging, OSC-TV was shown to outperform the analytic FDK approach and POCS-TV, an algorithm based on Gaussian statistics [27], in terms of image quality. Hence, it is conjectured that it may also benefit the 4D CBCT reconstruction. Its mathematical and physical bases have already been described in previous work, but it is still relevant to outline some of the theory supporting the algorithm. It has to be noted that this optimization assumes a mono-energetic photon beam and detector readings following a Poisson distribution.

The OSC step is performed first and attempts to maximize the Poisson log-likelihood L (μ) of the image estimate defined as follows:

$L (μ) = - \sum_{i} (d_{i} e^{- t_{i}} + Y_{i} t_{i}),$ (1) where μ refers to the voxels’ attenuation coefficient, i to a particular detector reading index, d_i to its incident photon count, Y_i to its detected photon count and t_i to the total attenuation along the ray path from the source to that detector. The dependency of L (μ) to μ is more straightforward when the mathematical definition of t_i is explicit:

$t_{i} = \sum_{j} l_{ij} μ_{j},$ (2) where l_ij represents the discrete length path in each voxel j. To maximize Equation 1, the convex algorithm [13] updates μ_j at each iteration n according to the following equation:

$μ_{j}^{(n + 1)} = μ_{j}^{(n)} + μ_{j}^{(n)} \frac{\sum_{i} l_{ij} [{\bar{y}}_{i}^{(n)} - Y_{i}]}{\sum_{i} l_{ij} t_{i}^{(n)} {\bar{y}}_{i}^{(n)}},$ (3) where ${\bar{y}}_{i}^{(n)} = d_{i} e^{- t_{i}^{(n)}}$ represents the expected photon count for a given detector and iteration. To accelerate the reconstruction, Kamphuis et al. [11] suggested to organize projections in subsets and to update the attenuation coefficients after each subset is considered. The more frequent update of μ does introduce a slight bias, but this bias is reduced by diminishing the number of subsets at the end of the reconstruction [2]. In the OSC-TV technique, the number of subsets is rather updated gradually. The number of subsets S for each iteration n is determined by Eq. 4 and depends of the initial number of subsets S_i, the final number of subsets and a reduction exponent p, with a typical value of p = 1/2. [16].∥ $| S |^{(n)} = round [\frac{| S |_{i} - | S |_{f}}{(n_{max} - 1)^{p}} (n_{max} - 1 - n)^{p} + | S |_{f}]$ (4) ∥Each subset contains more or less regularly spaced projection views. The subsets are selected in a way that maximizes the angular distance between the starting projection view of two consecutive subsets. Once the attenuation of each voxel is estimated by the OSC step, the TV step attempts to minimize ||μ||_TV, defined as the sum of the magnitudes of the discrete gradients of the image:

$| | μ | |_{TV} = \sum | \vec{\nabla} μ |,$ (5) TV regularization is performed iteratively by subtracting the local gradient for each voxel:

$μ_{j}^{(q + 1)} = μ_{j}^{(q)} - v_{j}^{(q)} \frac{c d_{A}}{\sqrt{\sum_{j} v_{j}^{2}}},$ (6) where q represents the TV iteration number, $v_{j} = \frac{δ | | μ | |_{TV}}{δ μ_{j}}$ the local value of the TV gradient for each voxel, $d_{A} \equiv \sqrt{\sum_{j} (μ_{j}^{(n)} - μ_{j}^{(n - 1)})^{2}}$ the norm of the difference between two OSC iterations, and c an empirical constant [27]. In the present work, 20 TV iterations were completed after each OSC iteration (with the exception of the first OSC iteration, for which d_A is not defined and therefore TV regularization is impracticable). A value of c = 0.2 was experimentally found to be suitable for 4D reconstructions, for both numerical phantom and clinical datasets.

The distance elements for the re-projection and backprojection were obtained via the Siddon’s method [26].

2.2 Adaptation to 4D reconstruction

The Amsterdam Shroud (AS) algorithm [37], as implemented in the Reconstruction Toolkit library (RTK) [20], was used to characterize the respiratory movement from a full cone beam CT acquisition. Though not as robust as the use of an external surrogate, the AS algorithm allows for 4D reconstruction without the need of additional apparatus for the patient. This signal was used for the identification of respiratory cycles and the following separation of the projections set into 8 respiratory subsets, each containing about 85 projections. In this work, the end expiration and end inspiration will respectively be defined as the 0% and 100% phases.

The tumor movement was assumed to be hysteresis-free and projections were sorted only according to the diaphragm relative displacement. Thus, projections acquired during inhalation and exhalation were grouped together given that they had the same diaphragm position. This type of separation could be source of position’s imprecisions for some patients, considering that tumor trajectory shows hysteresis for about 50% of the patients and that this hysteresis is superior to 2 mm for about 20% [24]. However, it is believed that such a phase separation method would allow for a better image quality for most patients. The reconstruction algorithm could also be applied to other types of respiratory subsets, but it was not investigated in the present study.

Three methods have been developed to adapt the OSC-TV algorithm to 4D reconstruction, each method having a different way of initializing the reconstructed image or, in mathematical terms, the attenuation coefficients of each voxel for the first iteration ( $μ_{j}^{(0)}$ in Eq. (3)). In the first method (4D OSC), a standard initialization of the volume grid is performed with a blank image. Each respiratory subset is therefore reconstructed completely separately from the other ones. In the second and third methods (3D OSC + 4D OSC and FDK + 4D OSC), the 4D reconstruction images are initialized from a prior 3D reconstruction, which can be performed by the OSC or FDK algorithm.

This prior reconstruction enables the use of information from the whole projections set, which could potentially reduce streaking artifacts in the final image. However, motion artifacts present in the initial 3D reconstruction could affect the final 4D image. We conjectured that those artifacts will be reduced by limiting the number of 3D iterations: while being blurrier, the initialization image is expected to be less prone to streaking artifacts after fewer iterations.

Figure 1 resumes the followed steps of the different methods for a 4D OSC-TV reconstruction.

Fig.1

Flowchart of the different proposed methods to adapt the OSC-TV algorithm to 4D reconstruction.

2.3 GPU implementation and hardware

The GPU’s parallel computing architecture is particularly suitable for massively parallel problems. Most of the different steps previously mentioned (raytracing, forward projection, back projection and TV regularization) can be separated in small tasks repeated multiple times (for example, for each voxel or each ray) and can therefore take advantage of a GPU implementation. The CUDA architecture (NVIDIA, Santa Clara, CA) was used to that end. Both reconstruction algorithms (FDK and OSC-TV) were implemented on GPU. The FDK implementation supplied by the RTK library was used [20].

In addition, the independent nature of each phase’s reconstruction facilitates the use of multiple GPUs. Considering the design constraints limiting the performance of a single GPU, multi-GPU computing clusters have recently emerged as leading supercomputing platforms, notably in the medical physics field [10]. Multi-GPU acceleration is particularly suitable for easily separable problems which require minimal inter-GPU communication. Indeed, once the optional prior 3D reconstruction and the respiratory signal are computed, respiratory subsets are reconstructed separately from each other, from different projections and a slightly different geometry. The 4D portion of the algorithm was therefore easy to implement on multiple GPUs, requiring no communication or synchronization between devices. It should be noted that the 3D portion of the algorithm was performed on a single GPU.

The reconstructions were performed on a node of 8 Tesla K20 GPUs (NVIDIA, Santa-Clara, CA) and 2 Intel Xeon processors (Intel, Santa Clara, CA). Tesla K20s are fitted with 2496 computing cores and have a global random access memory (RAM) size of 5 GiB for each GPU.

2.4 Projection data

The proposed algorithm was assessed both on the anthropomorphic numerical phantom XCAT [23] and on a series of clinical datasets.

The use of a numerical phantom was necessary to simulate a fully controllable acquisition, while the choice of this particular phantom was motivated by its accurate representation of human anatomy and respiratory movement. A total of 56 3D phantoms were generated to sample the respiratory movement, which had a period of 5 seconds (3 seconds of inspiration and 2 seconds of expiration), an amplitude of 1.2 cm in the anterior-posterior axis and of 2 cm in the longitudinal axis. The tested period and amplitude were respectively slightly longer and larger than most clinical data [22], but the ability to perform a valid reconstruction in such conditions would demonstrate a robustness to variable inter-patient’s movement patterns. A sphere of soft tissue of 1 cm diameter was also added in the right lung to simulate a lung tumor and to study the movement of small objects. Cardiac motion was not simulated.

The projections were obtained using the XCAT projector with parameters inspired from those of a Varian on-board imager (Palo Alto, CA) low-dose thorax scan. Attenuation of 70 keV photons was simulated for a half-fan detector of 397×298 mm², to obtain 672 projections on a fine grid of 1024×768 pixels. Two projection datasets of the numerical phantom were studied: one without noise and one with a level of noise matching the clinical projections described in the next paragraph. In order to reduce reconstruction time, the simulated projections were brought to 512×384 pixels via averaging of 2×2 pixel groups. The reconstruction grid consisted of 384×384×188 voxels of 1.2×1.2×1.5 mm³.

The clinical CBCT projection data of a patient, acquired by the on-board imager of Varian Clinac iX accelerators, were also used to investigate the performance of the different proposed methods. The selected scan had a visible diaphragm to facilitate the detection of the respiratory signal by the AS algorithm. Scans where the diaphragm is not visible at all times could be reconstructed with the proposed algorithm, but would necessitate other projections sorting methods, such as image intensity analysis [12] or the use of an external motion tracking system.

The clinical projections set contained 699 projections acquired in half-fan mode. The different reconstruction methods were tested on projections with full resolution (1024×768 pixels) and with reduced resolution (512×384 pixels). Due to GPU memory limitations, the 3D OSC-TV reconstruction from projections in full resolution was performed from half the projections (350 projections covering a full angular range). It was found that such a reconstruction demonstrated sufficient image quality, which was corroborated by a previous study [16].

2.5 Comparative evaluation

The 4D OSC-TV-based methods were compared to two well-known reconstruction algorithms in 4D CBCT: the MKB algorithm and the ASD-POCS algorithm. The MKB method is a 4D algorithm developed in the 1980s for the imaging of cardiac movement [17] and then adapted more recently for thoracic 4D CBCT [15]. An FDK reconstruction was first performed and used to obtain forward projections including motion. As suggested by Leng et al. [15], a median filter was then applied to the difference between initial projections and forward projected data to reduce streaking artifacts. Its kernel size was of 3×3 for 512×384 projections and of 5×5 for 1024×768 projections. This difference image was finally reconstructed and this reconstruction subtracted from the original reconstruction to obtain the final MKB image.

The ASD-POCS algorithm [28], on the other hand, is an iterative algorithm based on the minimization of the total variation of the image. It includes an algebraic reconstruction technique step with a positivity constraint and a reduction of the total variation by a gradient descent step. As proposed by Bergner et al. [4] and in a similar way to FDK + 4D OSC, a prior FDK reconstruction was used to initialize the respiratory reconstructions. For a fair comparison and a similar computational cost, the same number of iterations (projections-backprojections) was performed for the ASD-POCS reconstruction and for the 4D OSC-TV. The algorithm parameters for the numerical phantom with noiseless projections were based on the adaptation for 4D CBCT described by Bergner et al. [4], while those for the numerical phantom with noisy projections and for the clinical data were based on a clinical study performed by Shieh et al. [25].

The error of each reconstruction algorithm was estimated by the normalized root-mean-square deviation (NRMSD), defined as follows:

$NRMSD \equiv (\frac{1}{μ_{p, max} - μ_{p, min}}) \sqrt{\frac{\sum_{j} {(μ_{j, p} - μ_{j, r})}^{2}}{j_{max}}},$ (7) where μ_p denotes phantom voxel values and μ_r reconstructed voxel values. The normalization was performed with respect to the phantom voxel values extrema. NRMSD was computed in 3D for all the slices completely within the field of view (FOV). This metric may appear rather simple; nevertheless, it is of critical importance for prospective treatment techniques. In fact, 4D treatment planning could benefit from 4D CBCT [29, 33], where the target volume and the organs at risk should be segmented and registered to produce an accurate dynamic treatment plan. Therefore, the errors over the entire FOV should be considered in the reconstruction quality metric. In addition, this metric was computed for a 30 × 30 × 30 voxels region of interest (ROI) containing the lesion for all phases, to observe reconstruction error on the lesion movement trajectory.

For the clinical reconstructions, patient data was used, therefore, no reference measurements were available for the geometry or μ values. Patient images were analyzed from a qualitative standpoint, by observing noise, streaking artifacts and structure boundaries.

3 Results and Discussion

3.1 Simulated data

Projections of the moving phantom were first reconstructed by the two 3D algorithms, the OSC-TV iterative algorithm and the conventional FDK method (Fig. 2). For all reconstructed images (Figs. 2, 5 and 6), the same slice is shown with a μ range of [0, 0.3] cm^-1. In Fig. 2, both images display notable motion blurring and artifacts, which would compromise clinical use. These reconstructions are shown for comparison and to reaffirm the usefulness of a phase-correlated algorithm.

Fig.2

3D reconstruction of the XCAT phantom in movement using the FDK algorithm and the OSC-TV algorithm. Both reconstructions display notable motion blurring and artifacts. μ range of [0, 0.3] cm^-1 shown.

To optimize the FDK+4D OSC method, different backprojection filters were tested for the initializing 3D reconstruction. It was found that this parameter had little impact on the final 4D images, the subsequent OSC-TV iterations reducing efficiently the possible noise. For example, relative cutoff frequencies for the cosine-windowed ramp filter of 0.5, 0.75 and 1 were tested and yielded NRMSD values of 0.0467, 0.0471 and 0.0472 respectively, after 12 OSC-TV iterations for phase 0%. The cosine filter with a cutoff frequency of 1 was therefore selected for all the prior FDK reconstructions.

Figure 3 shows the progression of NRMSD as a function of the number of iterations of the different iterative algorithms. Phase 0% was selected to show the accuracy of reconstruction methods for a peak phase, which, as will be shown in Fig. 4, were more prone to artifacts. Remarkably, the type of prior influences the final result of the reconstruction. It is explained by the design features of OSC-TV, which uses alternating OSC and TV updates. The TV constraint is such that the final image estimate is in the vicinity, but not coinciding with the estimate which satisfies the OSC maximum likelihood. OSC-TV also uses convergence dampening via subset size reduction. Hence, the prior influences in which vicinity of the OSC optimum the estimation process will terminate.

Fig.3

NRMSD as a function of completed iterations for different iterative methods (3D OSC, 4D OSC, 3D OSC+4D OSC, FDK+4D OSC and ASD-POCS) for phase 0%. Results were obtained for the XCAT phantom using (a) noiseless and (b) noisy projection data. The four 4D methods yield a lower error than 3D OSC-TV, but the methods with a prior reconstruction yield the lowest error estimator and converged faster.

Fig.4

NRMSD as a function of the reconstructed phase for different methods (3D OSC, 4D OSC, 3D OSC+4D OSC, FDK+4D OSC, MKB and ASD-POCS). The NRMSD is measured for (a,c) the entire FOV and (b,d) a ROI of 30 × 30 × 30 voxels around the 1-cm lesion. With noiseless projection data (a,b), the four iterative methods yield a lower estimation error than MKB, while FDK+4D OSC and 3D OSC+4D OSC yield the lowest NRMSD. The same metrics were measured on the reconstruction of the XCAT with noise (c,d) and similar results were obtained.

Results of Fig. 3 were used to determine the number of iterations for each method. For the 3D OSC-TV reconstruction, NRMSD did not notably decrease after 2 iterations (NRMSD improvement inferior to 1%). Based on this result and after a few trials, it was decided to perform 2 iterations for the 3D part of the 3D OSC+4D OSC method. For the 4D OSC only method, a threshold of NRMSD improvement superior to 1% was selected as the best compromise between speed and accuracy. 12 iterations were therefore completed for this method. Considering the advantage provided by a prior image, the prior-initialized methods (3D OSC+4D OSC, FDK+4D OSC and ASD-POCS) were set to perform 10 respiratory-correlated iterations, whether it be performed by the OSC-TV or the ASD-POCS algorithms. In so doing, the different iterative methods perform a similar number of numerical operations.

Figure 4 presents the evolution of NRMSD with respect to the reconstructed phase for two different ROIs: the entire FOV (Fig. 4 (a)) and a ROI of 30 × 30 × 30 voxels around the 1-cm lesion (Fig. 4 (b)). This second ROI was selected to study the visualization of moving object more specifically. The different adaptations of the OSC-TV algorithm for 4D imaging can be compared between each other and with the MKB and ASD-POCS algorithms. It is clear that the OSC-TV methods yield a lower estimation error than the MKB method for all studied phases. The use of a prior 3D reconstruction also helps to slightly reduce the error estimate, especially for end-movement phases. ASD-POCS yields results similar to the other iterative algorithms for the full FOV, but represents the lesion ROI region less accurately. It is also worth noting that MKB yields a higher NRMSD than the 3D reconstruction performed by FDK for the full FOV, which can seem surprising. In fact, MKB is a method designed to visualize the movement, but also introduces artifacts in motionless regions and, in doing so, degrades the general image quality.

For all reconstruction methods, a lower NRMSD can be observed for the central phases. This decrease is attributable to the more even angular distribution of the projections used for those reconstructions. Indeed, projections identified as belonging to extremal movement phases, such as peak inhaled and peak exhaled, are bundled in one position of the respiratory cycle, while projections belonging to the middle movement phases appeared at two positions of the cycle, leading to a better angular coverage for middle phases.

Figures 3 and 4 also show the results obtained for the reconstruction of noisy projections. For the iterative methods, the NRMSD values after convergence are very similar whether or not the noise is simulated. The MKB and FDK reconstruction, however, were more affected by the noise. Noisy FDK initializations did not prevent the OSC 4D method from achieving the lowest error after convergence, as seen in Fig. 4 (c, d) The presence of noise was not visually identifiable on the reconstructed images themselves, so they were not included. Considering all these observations, the main source of error in the case of 4D CBCT reconstruction was the sparse and irregular angular sampling, while the imaging noise has much lower influence on image quality.

Figures 5 and 6 show, for phases 0% and 43%, the visual comparison between the different 4D reconstructions and the original phantom representing the central position of the studied phases. It can be observed that the use of a prior 3D reconstruction for initialization of the 4D reconstruction substantially increases image quality for an end phase such as phase 0%, while, for a central phase such as phase 43%, the three OSC-TV adaptations provide similar results.

Fig.5

XCAT phantom and its 4D reconstructions for respiratory phase 0% (end-expiration). The use of a prior image initialization considerably improved image quality. μ range of [0, 0.3] cm^-1 shown.

Fig.6

XCAT phantom and its 4D reconstructions for respiratory phase 43% (central respiration phase). Iterative methods provide comparable results, while MKB reconstruction shows both important artifacts and motion blurring. μ range of [0, 0.3] cm^-1 shown.

Qualitatively, 3D OSC+4D OSC and FDK+4D OSC reconstructions contain less undersampling artifacts in motionless regions (back and the exterior of the chest cavity) than the simple 4D OSC version, but some motion artifacts present in the 3D reconstruction are still visible, mostly in the sternum and heart regions. The ASD-POCS reconstruction, also initialized by a 3D image, is similar in that way to the 3D OSC+4D OSC and FDK+4D OSC reconstructions, but contains more noise and streaking in the uniform regions. On the other hand, the MKB reconstruction shows both important noise and motion blurring. If the displacement of a large object such as the diaphragm is noticeable, the visualization of the motion of small elements such as the tumor is much harder.

The temporal resolution and reduction of motion blurring were evaluated by the position detection of the simulated tumor and of the diaphragm. A separation in 8 respiratory subsets, assuming a hysteresis-free movement, allows for a theoretical average temporal resolution around 0.3 s (the half-period of 2.5 divided by the number of phases). Considering the amplitude of 1.2 cm in the anterior-posterior axis and of 2 cm in the longitudinal axis, a maximum theoretical residual movement about 0.15 cm (1.2 cm divided by the number of phases) in the anterior-posterior axis and about 0.25 cm (2 cm divided by the number of phases) in the longitudinal axis is expected. For both peak phases (phases 0% and 100%) and moving object (diaphragm in Fig. 8 and simulated tumor in Fig. 7), a profile was plotted, sampling a 3D volume plane with a line width of 3 pixels, along the anterior-posterior axis for the tumor and along the longitudinal axis for the diaphragm. As the tumor also moved in the longitudinal axis, its attenuation was studied in central slices of both peak phases. Position 0 approximately represents the center of the movement range.

Fig.7

Attenuation profiles for phase 0% (a) and 100% (b) as measured across the diaphragm in the longitudinal axis for different 3D volumes.

Fig.8

Attenuation profiles for phase 0% (top row) and 100% (bottom row), as measured where the 1 cm diameter soft tissue sphere is located during phase 0% (left column) and phase 100% (right column). Profiles are plotted in the anterior-posterior axis.

Figures 7 and 8 both show that motion blurring is largely reduced by the iterative methods. The study of the diaphragm movement demonstrates an accurate detection of its position and very low blurring around its edge. Attenuation coefficients match sufficiently well those of the XCAT phantom. In terms of tumor position, the three OSC-TV adaptations methods and ASD-POCS also provide similar results. The MKB reconstruction is the only reconstruction to show a clear contamination between phases: for the other methods, the simulated tumor is only visible on the expected slice. It should also be noted that FDK+4D OSC and 3D OSC+4D OSC demonstrate an attenuation coefficient which is more stable and closer to its actual value for both studied phases in the tumor region. Such attributes could facilitate a tumor delimitation process in a clinical context.

It has to be noted that the reconstructions tend to underestimate the tumor size in the posterior direction, of 0.12 cm for phase 0% and 0.24 cm for phase 100%. Such a slight offset can also be observed on the difference image presented in Figs. 5 and 6, which suggest a geometrical misalignment unrelated to the detection of movement. It could be explained by a possible small discrepancy between the continuous version of the XCAT phantom used to obtain projections and its voxelized version used to perform image quality analysis. However, for phase 100%, this could be exacerbated by a contamination from other respiratory phases.

The low tumor edge blurring and relatively precise tumor position detection suppose an adequate temporal resolution. However, detection of some of the real tumors could be complexified by their eventual proximity to artifacts. Since the initialization from a prior reconstruction reduces artifacts for end movement phases, detection of small objects’ motion could be facilitated by this initialization scheme.

From a more general perspective, it is understood that, by its simple movement (i.e., a respiratory movement of constant frequency and amplitude and no cardiac movement) and by the absence of noise in the data, the numeric phantom was a simplified model. Nevertheless, these conditions were already somewhat challenging and demonstrated the impact of the initialization image on convergence. The study of the 4D OSC-TV algorithm on clinical datasets is essential to demonstrate its interest.

3.2 Clinical data

A clinical acquisition was reconstructed successfully with the 4D OSC-TV methods. The 4D image allowed for the detection of a movement amplitude of 1.6 cm in the longitudinal axis, while no movement was observed in the anterior-posterior axis. Reconstructions of this patient with different algorithms are shown in Fig. 9.

Fig.9

Reconstructions of a clinical dataset (a) with the FDK algorithm (b) with various 4D algorithms for phases 0% and 100%. μ range of [0, 0.3] cm^-1 shown.

The four reconstruction methods show clearly the diaphragm movement in Fig. 9. A movement of smaller structures in the heart region is also visible on the OSC-TV-based reconstruction, but, in general, smaller structures seemed to be too regularized by the TV method for the clinical dataset. For example, the apparent lesion located in the lower region of the left lung (right side of the picture), suffers from loss of detail on small structures. It is still believed that, because of the high contrast between a lung tumor and the healthy surrounding tissues, a tumor’s movement could still be grossly visualized, but the loss of resolution is concerning. In addition to sparse angular sampling, it should be taken into consideration that the clinical dataset is quite ill-conditioned due to noise, scattered radiation, beam-hardening and gantry orbit irregularities typical to on-board CBCT. In consequence, iterative reconstruction mostly eliminates noise and high-frequency artifacts; in return, it suffers from high-frequency information loss, which is a limitation of the proposed technique consistent with other TV-based methods. Beyond the results presented here, smaller values of the TV regularization coefficient were also tested, but streaking artifacts due to irregular phase sampling were more prominent and the general image quality was not improved. The MKB clinical reconstruction is sharper, at the cost of added noise, which affects uniform structures as well as structure boundaries. For ASD-POCS, the streaking and noise levels were very sensitive to the value of regularization parameters and difficult to eliminate. These findings are consistent with the high-frequency noise visible in the 4D ASD-POCS reconstructions performed by Bergner et al. [4] and the persistent streaking reported by Matenine et al. [16] with respect to the 3D TV-regularized POCS method.

It appears clear that the image quality of reconstructions based on simulated data is noticeably better than the one of the clinical reconstruction. It is mainly attributed to the properties of simulated data, namely the absence of cardiac movement and the consistency of the respiratory movement, which all facilitated the task of the reconstruction algorithms. The TV regularization also performs best when filtering almost gradient-free regions such as those of the XCAT phantom. Furthermore, with respect to iterative reconstruction methods, the empirical parameters of the OSC-TV and ASD-POCS algorithms, such as the magnitude of the TV regularization, were straightforward to adjust for simulated data. For clinical data affected by noise and other movements, parameters became quite difficult to adjust and even there, the uncertainties which affect the acquisition were transferred to the reconstructed images.

The developed method was also validated on a reduced version of the clinical projection set, using a detector grid numerically binned to 512×384 pixels. The results are presented on Fig. 10 and are similar to the reconstruction from a full-resolution detector grid. It appears that for the geometry used in this work, the resolution of the detector is of relatively low importance when compared to image degradation due to sparse and irregularly spaced projection datasets.

Fig.10

Reconstructions of a clinical dataset with numerically reduced detector resolution using (a) the FDK algorithm and (b) different 4D methods for phases 0% and 100%. μ range of [0, 0.3] cm^-1 shown. The results are virtually identical to those obtained at full detector resolution.

3.3 Reconstruction times

To evaluate the execution time of the XCAT phantom, 5 repetitions of each reconstruction were performed to obtain a mean execution time and the corresponding standard deviation. For the 3D initialization images obtained using 3D OSC-TV, each iteration completed with the full dataset took 168.9±0.3 seconds. It is to be noted that 3D OSC-TV was performed on a single GPU. In contrast, the 4D reconstruction has two accelerating factors: the set of projections is about 8× smaller for each phase and the reconstruction of phases is performed concurrently on 8 GPUs. An iteration for the reconstruction of the the 8 phases concurrently on the 8 GPUs took 22.37±0.08 seconds on average. ASD-POCS iterations were slightly faster with an average computation time of 21.6±0.2 seconds per iteration. For its part, the FDK reconstruction took about 14±2 seconds. The total reconstruction times are shown in Table 1. Similar reconstruction times were obtained with the clinical dataset for a similar detector resolution of 512×384 pixels. The latter required between 5 % and 10 % more time than the XCAT phantom reconstruction. It is also worth noting that the MKB algorithm was not optimized and that most of its long computation time can be explained by the use of a median filter applied via CPU computations (10 min 8 sec).

Table 1
Reconstruction times of the XCAT phantom by the different 4D algorithms evaluated in this paper (for all 8 respiratory phases). Projection data: 672 views of 512×384 pixels; image grid: 384×384×188 voxels

Methods Computation times (s)

3D portion 4D portion Total

4D OSC - 298±11 298±11

FDK + 4D OSC 14±2 259±10 272±12

3D OSC + 4D OSC 337.9±0.7 248±9 586±9

MKB - ≈789 ≈789

ASD-POCS 14±2 216±2 230±4

Methods	Computation times (s)
4D OSC	-	298±11	298±11
FDK + 4D OSC	14±2	259±10	272±12
3D OSC + 4D OSC	337.9±0.7	248±9	586±9
MKB	-	≈789	≈789
ASD-POCS	14±2	216±2	230±4

It has to be noted that the 3D OSC+4D OSC adaptation is disadvantaged, in terms of computation time, by its use of a single GPU for the 3D OSC portion. Parallelization on multiple GPUs of the standard 3D OSC-TV code could reduce the computation time of 3D iterations up to a factor 8, making the total reconstruction time of 3D OSC + 4D OSC comparable to the 4D OSC method. However, FDK would remain faster than the iterative algorithm. Another advantage of the FDK initialization is the fact that, in a clinical context, the FDK reconstruction could begin during the scan acquisition. In such a case, the reconstruction time of the FDK + 4D OSC approach would only consist of its 4D portion.

Table 1 also shows that the ASD-POCS is slightly faster than the OSC-TV adaptations: for example, FDK+4D OSC take 15% more time to perform a 4D reconstruction than ASD-POCS. Both algorithms can therefore be interesting to use in an environment requiring a fast solution. However, it is important to remember that ASD-POCS has more parameters to adjust, which can complicate clinical software calibration. A different approach to acceleration was proposed by Zhong et al. with a mesh-based organ motion modeling, which in turn constrains a simultaneous motion estimation technique used in an iterative few-view 4D CBCT reconstruction algorithm [36]. Under these particular constraints, the algorithm is capable of operating on very sparse projection sets (about 20 views per phase), which contributes to the low reconstruction time of about 4 min for a 200×200×100 volume on a single NVIDIA Tesla K40c. In contrast, the method proposed here operates on about a hundred projections per phase, which requires a multi-GPU setup.

Computation times in the order of a few minutes are particularly interesting for clinical applications. Such reconstruction times are still too long to allow for the modification of a treatment plan with the patient on hold. However, a target movement model may be extracted from the 4D CBCT to characterize the treatment delivery and adjust the next treatment. For comparison purposes, the Table 2 shows various computation times and the corresponding projection and image grid sizes for a few algorithms of similar purpose. To our knowledge, even if some 4D reconstruction algorithms achieve acceptable results in less than 30 minutes [18 , 34] by parallelizing the computation on GPU, very few are faster than the 5 minutes mark [35]. It is evident that the implementation on multiple GPUs greatly reduced reconstruction times and such an approach could be used to accelerate other algorithms. However, parallelization on multiple GPUs would be more complex for algorithms requiring frequent communication between devices. The initialization by an analytical 3D reconstruction appears like a fast and simple way of increasing image quality by the use of knowledge from the full set of projections without compromising the ease of parallelization of each phase’s reconstruction. Without being real-time in the current form, the proposed solution can easily keep pace with the clinical patients flow and may be successfully used to visualize the patient’s anatomy.

Table 2

Computation times, as well as projection data and image grid sizes of comparable reconstruction algorithms

Authors	Computation time	Number of projections (size)	Number of reconstructed voxels
This work (numeric phantom)	4 min 32 sec	672 (512×384)	384×384×188
This work (patient data)	18 min	699 (1024×768)	384×384×188
Christoffersen et al. [6]	6 h	672 (512×384)	256×256×160
Tian et al. [30]	65 s	500	One axial slice
Gao et al. [9]	Less than 10 min	400 (256×64)	256×256×64
Zheng et al. [35]	20 s after the end of the acquisition	(504×376)	400×400×158
Yanet al. [34]	13 min	(512×384)	Not mentioned
Park et al. [19]	8 min	670 (512×512)	256×256×256
Mory et al. [18]	23 min	672 (512×384)	192×192×94

4 Conclusion

In the present study, the ability of the OSC-TV algorithm to perform 4D reconstructions on cone beam CT projection datasets obtained via the Amsterdam Shroud algorithm has been investigated. It was shown that the initialization reconstruction of each phase by a prior 3D reconstruction improved overall image quality. The 3D + 4D approaches were more robust, being able to reconstruct end movement phases with less artifacts in motionless regions, but showed slightly more streaking artifacts than the simple 4D approach. Both prior 3D initializations showed similar results, but the FDK initialization was significantly faster, with a total computation time of 272 s. Compared to the 4D ASD-POCS algorithm, the FDK+4D OSC method was shown to reconstruct a moving object more accurately, at the cost of 15% of added reconstruction time.

Compared to many 4D algorithms and their numerical benchmarks available in the literature, FDK + 4D OSC-TV compares favourably in terms of computation time, thanks to the implementation on multiple GPUs. One of the purposes of using CBCT imaging for IGRT is to accurately and precisely identify the tumor position of the day as well as its trajectory during breathing. In this context, it is important to reduce motion blurring as much as possible. 4D CBCT can be used for that purpose, but only if image reconstruction can be performed rapidly. Our work shows that this is indeed feasible to achieve without large compromises on image quality nor spatial resolution. In addition, results were obtained from a CBCT acquisition with a standard rotation speed and number of projections, which would also facilitate a clinical implementation. It implies that the algorithm could be applied a posteriori on pre-existing scans if necessary.

Further work is focused on obtaining a better compromise between image quality and speed. The number of iterations, the importance of regularization and the reconstruction grid size could be optimized to correspond appropriately to clinical needs. For example, due to a very small movement in a certain direction, reconstruction’s voxels could be expanded in this direction to reduce computation time. The use of a 4D TV regularization could also improve image quality with a low computational cost, if adapted adequately to multiple GPUs.

Footnotes

Acknowledgments

The authors acknowledge partial support by the CREATE Medical Physics Research Training Network grant of the Natural Sciences and Engineering Research Council of Canada (Grant number: 432290).

Computations were made on the supercomputer Helios from Université Laval, managed by Calcul Québec and Compute Canada. The operation of this supercomputer is funded by the Canada Foundation for Innovation (CFI), the Ministère de l’Économie, de la science et de l’innovation du Québec (MESI) and the Fonds de recherche du Québec - Nature et technologies (FRQ-NT).

References

Balter

, Haken

R.K.T.

, Lawrence

T.S.

et al., Uncertainties in CT-based radiation therapy treatment planning associated with patient breathing, International Journal of Radiation Oncology Biology Physics 36 (1996), 167–174.

Beekman

F.J.

and Kamphuis

, Ordered subset reconstruction for x-ray CT, Physics in Medicine and Biology 46 (2001), 1835–1844.

Bergner

, Berkus

, Oelhafen

et al., Autoadaptive phase-correlated (AAPC) reconstruction for 4D CBCT, Medical Physics 36 (2009), 5695–5706.

Bergner

, Berkus

, Oelhafen

et al., An investigation of 4D cone-beam CT algorithms for slowly rotating scanners, Medical Physics 37 (2010), 5044–5053.

Chen

G.H.

, Tang

and Leng

, Prior image constrained compressed sensing (PICCS): A method to accurately reconstruct dynamic CT images from highly undersampled projection data sets, Medical Physics 35 (2008), 660–663.

Christoffersen

C.P.V.

, Hansen

, Poulsen

et al., Registration-Based Reconstruction of Four-Dimensional Cone Beam Computed Tomography, IEEE Transactions on Medical Imaging 32 (2013), 2064–2077.

Feldkamp

L.A.

, Davis

L.C.

and Kress

J.W.

, Practical cone-beam algorithm, J Opt Soc Amer A (1984), 612–619.

Gang

G.J.

, Stayman

J.W.

, Ehtiati

et al., Task-driven image acquisition and reconstruction in cone-beam CT, Physics in Medicine and Biology 60 (2015), 3129.

Gao

, Guntupalli

R.K.

, Lin

et al., 4D cone beam CT via spatiotem-poral tensor framelet, Medical Physics 39 (2012), 6943–6946.

10.

Jia

and Jiang

, eds., Graphics Processing Unit-Based High Performance Computing in Radiation Therapy, CRC Press;, 2015.

11.

Kamphuis

and Beekman

F.J.

, Accelerated iterative transmission CT reconstruction using an ordered subsets convex algorithm, IEEE Transactions on Medical Imaging 17 (1998), 1101–1105.

12.

Kavanagh

, Evans

P.M.

, Hansen

V.N.

et al., Obtaining breathing patterns from any sequential thoracic x-ray image set, Physics in Medicine and Biology 54 (2009), 4879–4888.

13.

Lange

and Fessler

J.A.

, Globally convergent algorithms for maximum a posteriori transmission tomography, IEEE Transactions on Image Processing 4 (1995), 1430–1438.

14.

Leng

, Tang

, Zambelli

et al., High temporal resolution and streak-free four-dimensional cone-beam computed tomography, Physics in Medicine and Biology 53 (2008), 5653–5673.

15.

Leng

, Zambelli

, Tolakanahalli

et al., Streaking artifacts reduction in four-dimensional cone-beam computed tomography, Medical Physics 35 (2008), 4649–4659.

16.

Matenine

, Goussard

and Després

, GPU-accelerated regularized iterative reconstruction for few-view cone beam CT, Medical Physics 42 (2015), 1505–1517.

17.

McKinnon

G.C.

and Bates

R.H.T.

, Towards imaging the beating heart usefully with a conventional CT scanner SCANNER, IEEE Transactions on Biomedical Engineering 28 (1981), 123–127.

18.

Mory

, Auvray

, Zhang

et al., Cardiac C-arm computed tomography using a 3D+time ROI reconstruction method with spatial and temporal regularization, Medical Physics 41 (2014), 021903.

19.

Park

J.C.

, Zhang

, Chen

et al., Common-mask guided image reconstruction (c-MGIR) for enhanced 4D cone-beam computed tomography, Physics in Medicine and Biology 60 (2015), 9157.

20.

Rit

, Oliva

M.V.

, Brousmiche

et al., The Reconstruction Toolkit (RTK), an open-source cone-beam CT reconstruction toolkit based on the Insight Toolkit (ITK), Journal of Physics: Conference Series 489 (2014), 012079.

21.

Ritschl

, Sawall

, Knaup

et al., Iterative 4D cardiac micro-CT image reconstruction using an adaptive spatio-temporal sparsity prior, Physics in Medicine and Biology 57 (2012), 1517.

22.

Schmidt

M.L.

, Poulsen

, Toftegaard

et al., Clinical use of iterative 4D-cone beam computed tomography reconstructions to investigate respiratory tumor motion in lung cancer patients, Acta Oncologica 53 (2014), 1107–1113.

23.

Segars

W.P.

, Mahesh

, Beck

T.J.

et al., Realistic CT simulation using the 4D XCAT phantom, Medical Physics 35 (2008), 3800–3808.

24.

Seppenwoolde

, Shirato

, Kitamura

et al., Precise and real-time measurement of 3D tumor motion in lung due to breathing and heartbeat, measured during radiotherapy, International Journal of Radiation Oncology Biology Physics 53 (2002), 822–834.

25.

Shieh

, Kipritidis

, O’Brien

R.T.

et al., Image quality in thoracic 4D cone-beam CT: A sensitivity analysis of respiratory signal, binning method, reconstruction algorithm, and projection angular spacing, Medical Physics 41 (2014), 041912.

26.

Siddon

R.L.

, Fast calculation of the exact radiological path for a threedimensional CT array, Medical Physics 12 (1985), 252–255.

27.

Sidky

E.Y.

, Kao

C.M.

and Pan

, Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT, X-ray Sci Tech 14 (2006), 119–139.

28.

Sidky

E.Y.

and Pan

, Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization, Physics in Medicine and Biology 53 (2008), 4777.

29.

Sweeney

R.A.

, Seubert

, Stark

et al., Accuracy and inter-observer variability of 3D versus 4D cone-beam CT based image-guidance in SBRT for lung tumors, Radiation Oncology 7 (2012), 81.

30.

Tian

, Jia

, Dong

et al., Low-dose 4DCT reconstruction via temporal nonlocal means, Medical Physics 38 (2011), 1359–1365.

31.

Wang

and Gu

, High-quality four-dimensional cone-beam CT by deforming prior images, Physics in Medicine and Biology 58 (2013), 231–246.

32.

Wang

and Gu

, Simultaneous motion estimation and image reconstruction (SMEIR) for 4D cone-beam CT, Medical Physics 40 (2013), 101912.

33.

Wijesooriya

, Weiss

, Dill

et al., Quantifying the accuracy of automated structure segmentation in 4D CT images using a deformable image registration algorithm, Medical Physics 35 (2008), 1251–1260.

34.

Yan

, Zhen

, Folkerts

et al., A hybrid reconstruction algorithm for fast and accurate 4D cone-beam CT imaging, Medical Physics 41 (2014), 071903.

35.

Zheng

, Sun

, Pavkovich

, et al., Fast 4D Cone-Beam Reconstruction Using the McKinnon-Bates Algorithm with Truncation Correction and Non Linear Filtering, in: SIE Medical Imaging, International Society for Optics and Photonics 7961 (2011), 79612U.

36.

Zhong

, Gu

, Mao

et al., 4D cone-beam CT reconstruction using multi-organ meshes for sliding motion modeling, Physics in Medicine and Biology 61 (2016), 996.

37.

Zijp

L.J.

, Sonke

J.J.

and van Herk

, Extraction of the respiratory signal from sequential thorax Cone-Beam X-ray images, in: International Conference on the Use of Computers in Radiation Therapy, May 10–13, (2004), Seoul, Korea, 507–509.

Methods	Computation times (s)
	3D portion	4D portion	Total
4D OSC	-	298±11	298±11
FDK + 4D OSC	14±2	259±10	272±12
3D OSC + 4D OSC	337.9±0.7	248±9	586±9
MKB	-	≈789	≈789
ASD-POCS	14±2	216±2	230±4