1. Introduction
Four dimensional cone-beam CT (CBCT) is required for many clinical tasks in radiotherapy, such as image guidance, 4D dose reconstruction, marker-less tumor tracking, and accurate volume assessment of moving tumors. It has also been demonstrated to reduce setup errors when compared with normal 3D CBCT [
[1]- Sweeney R.A.
- Seubert B.
- Stark S.
- et al.
Accuracy and inter-observer variability of 3D versus 4D cone-beam CT based image-guidance in SBRT for lung tumors.
]. However, long scan times and poor image quality remain a challenge for its widespread clinical use. Additionally, dedicated 4D protocols will typically have poorer 3D image quality, due to reduced tube current [
[2]- Thengumpallil S.
- Smith K.
- Monnin P.
- Bourhis J.
- Bochud F.
- Moeckli R.
Difference in performance between 3D and 4D CBCT for lung imaging: a dose and image quality analysis.
]. Iterative cone-beam CT reconstruction has been shown to allow for reduced imaging dose and improved image quality [
[3]Iterative image reconstruction for CBCT using edge-preserving prior.
]. For 4D CBCT, it can dramatically reduce scan times [
[4]- Qi Z.
- Chen G.H.
- Rit S.
- et al.
Performance studies of four-dimensional cone beam computed tomography.
]. On the other hand, reconstruction times are quite long ranging from 5 min to several hours in literature [
[5]- Christoffersen C.P.V.
- Hansen D.
- Poulsen P.
- Sørensen T.S.
Registration-based reconstruction of four-dimensional cone beam computed tomography.
]. Additionally, many algorithms require manual steps, such as segmenting the lung, which pose a challenge in a clinical setting. Recently several groups have demonstrated impressive results by using registration based methods for 4D CBCT [
5- Christoffersen C.P.V.
- Hansen D.
- Poulsen P.
- Sørensen T.S.
Registration-based reconstruction of four-dimensional cone beam computed tomography.
,
6- Harris W.
- Zhang Y.
- Yin F.F.
- Ren L.
Estimating 4D-CBCT from prior information and extremely limited angle projections using structural PCA and weighted free-form deformation for lung radiotherapy.
,
7- Zhang Y.
- Yin F.F.
- Segars W.P.
- Ren L.
A technique for estimating 4D-CBCT using prior knowledge and limited-angle projections.
,
8- Mory C.
- Janssens G.
- Rit S.
Motion-aware temporal regularization for improved 4D cone-beam computed tomography.
,
9- Liu J.
- Zhang X.
- Zhang X.
- et al.
5D respiratory motion model based image reconstruction algorithm for 4D cone-beam computed tomography.
,
10High-quality four-dimensional cone-beam CT by deforming prior images.
,
11Simultaneous motion estimation and image reconstruction (SMEIR) for 4D cone-beam CT.
,
12- Zhong Z.
- Gu X.
- Mao W.
- Wang J.
4D cone-beam CT reconstruction using multi-organ meshes for sliding motion modeling.
]. All of these methods are however based on custom scan sequences taking anywhere between 2 and 8 min. Reconstruction times range between 4 min and 15 h, and combined the fastest method takes about 9 min from scan start to reconstructed image [
[12]- Zhong Z.
- Gu X.
- Mao W.
- Wang J.
4D cone-beam CT reconstruction using multi-organ meshes for sliding motion modeling.
]. This makes their routine use impractical in a standard clinical workflow. Ideally, 4D reconstructions would be available within a few minutes, based on a standard 60 s 3D acquisition. In this work we present such an algorithm for 4D CBCT reconstruction, which runs in less than 60 s on a modern Graphics Processing Unit (GPU), and a new regularizer for improved temporal resolution. These were compared with state of the art iterative methods. For all methods regularization parameters were found automatically, via surrogate optimization.
2. Materials and methods
Iterative CBCT reconstruction is often treated as a regularized least squares problem
where
is the projection transform,
is the (4D) image to be reconstructed and
are the log-transformed projections.
is either the identity transform or (as in this work) a diagonal matrix containing weights [
[13]Optimal short scan convolution reconstruction for fan beam CT.
] to compensate for the offset detectors. Finally,
is the (nonlinear) regularization function, such as total variation, given by
(2)
in the standard 3D isotropic version. Here
t is the phase index, indicating that the total TV is the sum of the TV for the individual phases. We investigated three different choices for regularization.
2.1 Sparse Frequency Regularization (SFR)
The following regulariser is proposed:
(3)
denotes a box filter, downsampling the image by a factor of two in the spatial dimensions.
is the Fourier transform along the temporal dimension.
,
and
are unitless weights, controlling the strength of each regularisation parameter. The total variation term ensures that the images are sparse spatially, while the Fourier term enforces that the change in each pixel can be represented with a few Fourier coefficients. The use of total variation on multiple scales was previously proposed by Huang et al. [
[14]Huang Y, Taubmann O, Huang X, Haase V, Lauritsch G, Maier A. A new scale space total variation algorithm for limited angle tomography. In: CT-Meeting 2016 proceedings (the 4th international meeting on image formation in X-ray computed tomography); 2016. p. 149–52.
] where a Gaussian filter was used as the scaling operator. Conceptually, the use of TV on several resolution scales is also similar to what is done for Spatiotemporal framelet (STF) [
[15]- Gao H.
- Li R.
- Lin Y.
- Xing L.
4D cone beam CT via spatiotemporal tensor framelet.
], but only the first-order differential was included, in the interest of retaining clinical reconstruction times.
2.2 4D total variation
Several previous publications have demonstrated the effective use of 4D TV for CT reconstruction [
5- Christoffersen C.P.V.
- Hansen D.
- Poulsen P.
- Sørensen T.S.
Registration-based reconstruction of four-dimensional cone beam computed tomography.
,
8- Mory C.
- Janssens G.
- Rit S.
Motion-aware temporal regularization for improved 4D cone-beam computed tomography.
]. In this work 4D TV with isotropic TV in the spatial domain (Eq. (
2)), and anisotropic in the temporal dimension was implemented. Anisotropic TV was used as it significantly reduces the memory requirements of the algorithm, as each 3D volume can be updated separately. The regularizer is then given by
(4)
where
is the temporal total variation, given by
(5)
Again,
and
are unitless weights, controlling strength of the regularisation.
2.3 Prior image constrained compressed sensing (PICCS)
PICCS [
[16]- Chen G.H.
- Tang J.
- Leng S.
Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets.
] is an extension of standard 3D TV. First, a 3D image
is reconstructed using the standard Feldkamp-Davis-Kress (FDK) [
[17]- Feldkamp L.A.
- Davis L.C.
- Kress J.W.
Practical cone-beam algorithm.
] method. This is then used as a prior, and the regularization term becomes
(6)
where
is the 3D image corresponding to phase
t of
. The motivation behind PICCS is that most of the image is stationary, so the difference term
should mainly contain zero. While the original article used a constrained approach Lauzier et al. [
[18]- Lauzier P.T.
- Tang J.
- Chen G.H.
Prior image constrained compressed sensing: implementation and performance evaluation.
] established that the method also worked well in an unconstrained setting. Again,
and
are unitless weights controlling strength of the regularisation. For this method, the prior image was also used as a starting guess, rather than the standard zero filled image.
2.4 McKinnon-Bates
The McKinnon-Bates (MKB) algorithm [
[19]- Mc Kinnon G.C.
- Bates R.H.
Towards imaging the beating heart usefully with a conventional CT scanner.
] is a non-iterative algorithm commonly used for 4D CBCT, and can be considered clinical standard. In the first step of the algorithm a standard 3D image is reconstructed using the FDK algorithm. For each time phase a correction to the 3D image is calculated to produce the 4D image. The correction is made by forwards projection the 3D image and subtracting the original projections to create error projections. The FDK algorithm is then used on the error projections to create the correction.
2.5 RTK ROOSTER
The RTK framework [
[20]- Rit S.
- Vila Oliva M.
- Brousmiche S.
- Labarbe R.
- Sarrut D.
- Sharp G.C.
The Reconstruction Toolkit (RTK), an open-source cone-beam CT reconstruction toolkit based on the Insight Toolkit (ITK).
] is another open-source framework for CBCT reconstruction, which is widely used in the field. We compare with the state of the art 4D method ROOSTER [
[21]Mory C, Jacques L. A modified 4D ROOSTER method using the Chambolle–Pock algorithm. In: Proc. 3rd intl. conf. on image formation in X-ray CT; 2014. p. 191–3.
] implemented in this framework. We do not use the motion-aware version [
[8]- Mory C.
- Janssens G.
- Rit S.
Motion-aware temporal regularization for improved 4D cone-beam computed tomography.
], as we do not have the deformable vector fields required for this method. The regularization parameters used were
and
, which were found manually. We used 100 iterations as it resulted in the best quality of images. We also included 10 iterations as in the original paper, for comparison. All possible GPU acceleration was enabled, to ensure the fastest runtime.
2.6 Patient scans
Nineteen consecutive lung cancer patients, treated with Stereotactic Body Radiation Therapy (SBRT) were scanned using the onboard imager on a Varian Trilogy linear accelerator as part of their treatment. For further details on these patients, see Schmidt et al. [
[22]- Schmidt M.L.
- Poulsen P.R.
- Toftegaard J.
- Hoffmann L.
- Hansen D.
- Sørensen T.S.
Clinical use of iterative 4D-cone beam computed tomography reconstructions to investigate respiratory tumor motion in Lung Cancer patients.
]. The scans consisted of 663 projections in
over 60 s, using an offset detector array for extended field of view. The low-dose thorax protocol was used, using a voltage of 110 kV, a current of 20 mA and a pulse-length of 20 ms.
2.7 Extraction of breathing curves
In order to extract the respiratory phases from the projections, each projection was registered to the next using the Demons deformable image registration algorithm [
[23]Image matching as a diffusion process: an analogy with Maxwell’s demons.
]. The mean deformation in the craniocaudal direction was then integrated to provide a respiratory signal. This was filtered to remove the low frequency contribution from the rotation of the gantry and the high frequency contribution from cardiac motion. Finally, the inhale peaks of the extracted breathing curve was used to sort the projections into 10 respiratory phases using phase binning. This method was chosen rather than the standard Amsterdam Shroud [
[24]Zijp L, Sonke JJ, van Herk M. Extraction of the respiratory signal from sequential thorax cone-beam X-ray images. In: International conference on the use of computers in radiation therapy; 2004. p. 507–9.
], as the diaphragm was not consistently visible in all projections causing that method to fail.
2.8 Virtual XCAT phantom
To validate our approach and tune the regularisation parameters the female XCAT 4D phantom [
[25]- Segars W.P.
- Sturgeon G.
- Mendonca S.
- Grimes J.
- Tsui B.M.W.
4D XCAT phantom for multimodality imaging research.
] was used. A 1 cm-diameter spherical lesion was added to the bottom half of the right lung. The default bone thickness in XCAT is considerably greater than what was seen in our patient cohort. Bone thickness was therefore reduced by 50% compared to the XCAT default. The setup was designed to match clinical reality as closely as possible. Distance between source and detector was 150 cm and distance from source to isocenter was 100 cm. The detector was
with an image resolution of
. The detector was offset with 144.97 mm, equal to what is used in the clinical scans. 620 projections were obtained over a
arc. Scan duration was 60 s, with a regular breathing period of 4 s, similar to our patients. The projections were binned into 10 respiratory phases.
Ground truth images were created with the XCAT program. One motion free image was created for each projection time point. The images in the same respiratory bin were then averaged to provide a ‘best case’ ground truth. The voxel size was which was rebinned to .
2.9 Implementation
To solve Eq. (
1), an accelerated ordered subset (OS) algorithm was used [
[26]- Hansen D.C.
- Sangild Sørensen T.
- Rit S.
Fast reconstruction of low dose proton CT by sinogram interpolation.
] (see Appendix). This is a modification of the algorithm proposed by Kim et al. [
[27]- Kim D.
- Ramani S.
- Fessler J.A.
Combining ordered subsets and momentum for accelerated X-ray CT image reconstruction.
], adapted to allow for isotropic total variation by combining it with the work of Arrow et al. [
[28]- Arrow K.J.
- Hurwicz L.
- Uzawa H.
Studies in linear and non-linear programming.
]. The approach is similar to the one recently presented by Xu et al. [
[29]- Xu Q.
- Yang D.
- Tan J.
- Sawatzky A.
- Anastasio M.A.
Accelerated fast iterative shrinkage thresholding algorithms for sparsity-regularized cone-beam CT image reconstruction.
] for 3D CBCT, but was developed independently. An important advantage of this algorithm is that memory use can be kept low, allowing for the entire reconstruction to be done exclusively on the GPU without transferring back and forth between GPU and system memory.
2.10 Selection of subsets
While OS algorithms are widely used for CT and 3D CBCT reconstruction, they have not previously been used in 4D CBCT, presumably due to the problem of convergence. OS algorithms split the data into a
N subsets, and then perform a step of an iterative algorithm using only this subset. As the cost of an iteration is usually proportional to the size of the data, this provides a theoretical
N times speed up. Recently, ordered subsets have been combined with Nesterov’s method [
27- Kim D.
- Ramani S.
- Fessler J.A.
Combining ordered subsets and momentum for accelerated X-ray CT image reconstruction.
,
30Introductory lectures on convex optimization: a basic course.
] for further speed up. The number of subsets used does however greatly affect the stability of the algorithm. As 4D CBCT is already tenfold undersampled (due to the 10 breathing phases), this makes the exact strategy for selection of the subsets quite important. In this work, projections were put into the subsets in a straightforward round-robin approach, based on projection angle. As projections from the same phase come in “bunches” of 3 or 4, this spreads the projections in each bunch into different subsets. A relatively small number of subsets was employed (
6) to ensure algorithmic stability. Further increasing the number of subsets would only yield a small improvement in performance, as the regularization step becomes the performance bottleneck.
2.11 Forward and backward projection
The GPU implementation of the basic forward and backward projection (
and
respectively) is based on our previously published work [
[5]- Christoffersen C.P.V.
- Hansen D.
- Poulsen P.
- Sørensen T.S.
Registration-based reconstruction of four-dimensional cone beam computed tomography.
], although significant optimizations were performed. For the forwards projection, each image phase was loaded separately into a 3D GPU texture. We used one GPU thread per projection pixel. Raytracing was then performed by evaluating the line from the X-ray source to the detector pixel at
n equally spaced points through the volume, using trilinear interpolation. For the backprojection, the projections belonging to each phase were loaded into a stacked 2D texture. A GPU thread was used per image voxel and backprojection was performed by finding the corresponding point on each projection and sampling with bilinear interpolation. This approach is similar to what is used in other GPU reconstruction frameworks such as RTK [
[20]- Rit S.
- Vila Oliva M.
- Brousmiche S.
- Labarbe R.
- Sarrut D.
- Sharp G.C.
The Reconstruction Toolkit (RTK), an open-source cone-beam CT reconstruction toolkit based on the Insight Toolkit (ITK).
].
2.12 Parameters and evaluation
All images were reconstructed at a resolution of
. 10 full iterations were chosen, based on previous experience with the same algorithm for other applications [
[26]- Hansen D.C.
- Sangild Sørensen T.
- Rit S.
Fast reconstruction of low dose proton CT by sinogram interpolation.
]. To evaluate image quality on the XCAT phantom, structural similarity (SSIM) was used [
[31]- Wang Z.
- Bovik A.
- Sheikh H.
- Simoncelli E.
Image quality assessment: from error visibility to structural similarity.
]. SSIM is a metric, ranging from zero to one, which evaluates the similarity of two images. We used SSIM in 3D with
(1.5 pixels), and report the minimum value over all phases. As background artifacts were not of interest in this study, the background is removed before evaluation, based on the ground truth. Based on this, for all algorithms, the parameters (
and the respective
’) were found using 500 iterations of surrogate optimization [
[32]Mixture surrogate models based on Dempster-Shafer theory for global optimization problems.
] using pySOT [
[33]Eriksson D, Bindel D, Shoemaker C. Surrogate Optimization Toolbox (pySOT) [Internet]. Cited [2018 Mar 2]. Available from: github.com/dme65/pySOT.
]. The found parameters were then used for the patient scans.
4. Discussion
In this study we have demonstrated the possibility of using 4D iterative reconstruction within one minute. On phantom data the fast iterative methods outperformed non-iterative methods significantly, and matched that of iterative algorithms more than 20 times slower. When using the SSIM optimized parameters for patient scans PICCS generally provided the sharpest images, but also the largest degree of motion artefacts, including streaks. This is particularly pronounced in patient 19, who also had the slowest breathing. The SFR method consistently provides sharper images than TV across the three patients, but without increased streaking and motion blur. That being said, the differences are generally modest, and all three methods provide similar image quality.
The ROOSTER method provided similar image quality to the iterative methods used in this work. This is unsurprising, as the fundamental mathematical problem solved in the ROOSTER method is the same that is used in the TV4D method. It is however too slow to be used as part of the daily CBCT scan prior to treatment, which is the main limitation of that particular method. The long reconstruction time also means that automatic parameter optimization is not feasible, and regularisation parameters thus have to be found manually.
For radiotherapy, it is important that the tumour position is unbiased. While verifying this is difficult without implanted markers, it should be noted that PICCS only depends on the 3D FDK image, and the projections in the corresponding phase. Artefacts would thus present themselves as blurring, rather than a position bias. For the other methods, they show no bias when compared with the PICCS methods.
Only a few papers have investigated 4D CBCT from standard 3D acquisitions, presumably due to the challenges in retaining good image quality. We have demonstrated that images which could be clinically useful can indeed be obtained, which would be adequate for patient positioning and evaluation of tumour size.
In patients, irregular breathing and general movement can cause image artefacts. Even though this was not included in the digital phantom, the parameters found through the automatic search translated well to the patient scans, even when the patient breathing cycle did not match the sinosoidal 4 s curve used for the phantom.
It is an open question if the same parameters would perform as well when used with different scan protocols or at different sites, however. For longer acquisition times, it is likely that weaker regularization may improve results, as the base image quality is better. Investigating this is left as an exercise for the reader.
It is worth noting that patient 19 was reported as a failure case in a previous work on iterative 4D CBCT [
[22]- Schmidt M.L.
- Poulsen P.R.
- Toftegaard J.
- Hoffmann L.
- Hansen D.
- Sørensen T.S.
Clinical use of iterative 4D-cone beam computed tomography reconstructions to investigate respiratory tumor motion in Lung Cancer patients.
], where the tumor could not be distinguished. That the tumor is clearly visible in this work, for both SFR and TV, which should mainly be ascribed to the use of automatic parameter tuning, which was enabled by the fast reconstruction time.
Many papers have investigated 4D CBCT, but most rely on longer acquisition times, ranging from 2 to 9 min in duration [
8- Mory C.
- Janssens G.
- Rit S.
Motion-aware temporal regularization for improved 4D cone-beam computed tomography.
,
11Simultaneous motion estimation and image reconstruction (SMEIR) for 4D cone-beam CT.
,
34- Bergner F.
- Berkus T.
- Oelhafen M.
- et al.
An investigation of 4D cone-beam CT algorithms for slowly rotating scanners.
,
35- Vergalasova I.
- Cai J.
- Yin F.F.
A novel technique for markerless, self-sorted 4D-CBCT: feasibility study.
,
36- Shieh C.C.
- Kipritidis J.
- O’Brien R.T.
- Kuncic Z.
- Keall P.J.
Image quality in thoracic 4D cone-beam CT: a sensitivity analysis of respiratory signal, binning method, reconstruction algorithm, and projection angular spacing.
].
A few other works have investigated 4D CBCT from 60 s acquisitions such as Yan et al. [
[37]- Yan H.
- Zhen X.
- Folkerts M.
- et al.
A hybrid reconstruction algorithm for fast and accurate 4D cone-beam CT imaging.
] and Zhong et al. [
[12]- Zhong Z.
- Gu X.
- Mao W.
- Wang J.
4D cone-beam CT reconstruction using multi-organ meshes for sliding motion modeling.
]. These works suffer from the problem that they “simulate” a 60 s acquisition by subsampling a longer (4–6 min) acquisitions. This results in an acquisition where the breathing phases are more evenly spread out between the projections than in a true 60 s acquisition – resulting in a higher image quality. Qi et al. [
[4]- Qi Z.
- Chen G.H.
- Rit S.
- et al.
Performance studies of four-dimensional cone beam computed tomography.
] and Gao et al. [
[15]- Gao H.
- Li R.
- Lin Y.
- Xing L.
4D cone beam CT via spatiotemporal tensor framelet.
] studied the performance of various algorithms for 4D CBCT on 60 s and 30 s data respectively, but in both cases only on digital phantoms. To our knowledge, the only studies on 60 s data was then from our own group, where Christoffersen et al. [
[5]- Christoffersen C.P.V.
- Hansen D.
- Poulsen P.
- Sørensen T.S.
Registration-based reconstruction of four-dimensional cone beam computed tomography.
] investigated a new registration based algorithm and Schmidt et al. [
[22]- Schmidt M.L.
- Poulsen P.R.
- Toftegaard J.
- Hoffmann L.
- Hansen D.
- Sørensen T.S.
Clinical use of iterative 4D-cone beam computed tomography reconstructions to investigate respiratory tumor motion in Lung Cancer patients.
] evaluated it on a larger clinical cohort. In comparison, this work achieves a subjectively better image quality and reconstruction times which are at least two orders of magnitude faster.
In conclusion, we have presented a fast method for iterative 4D CBCT running in less than 60 s, and used it to compare two different iterative methods with one new, based on Fourier analysis. We have demonstrated that regularization parameters found via automatic optimization can translate well into clinical cases, and the parameters are robust over a larger cohort of patients. The next stop would be to perform clinical implementation research to introduce this into clinical practice.
Article info
Publication history
Published online: March 08, 2018
Accepted:
February 15,
2018
Received in revised form:
February 9,
2018
Received:
August 10,
2017
Copyright
© 2018 The Authors. Published by Elsevier B.V. on behalf of European Society of Radiotherapy & Oncology.