Advertisement

Fast 4D cone-beam CT from 60 s acquisitions

Open AccessPublished:March 08, 2018DOI:https://doi.org/10.1016/j.phro.2018.02.004

      Abstract

      Background & purpose

      Four dimensional Cone beam CT (CBCT) has many potential benefits for radiotherapy but suffers from poor image quality, long acquisition times, and/or long reconstruction times. In this work we present a fast iterative reconstruction algorithm for 4D reconstruction of fast acquisition cone beam CT, as well as a new method for temporal regularization and compare to state of the art methods for 4D CBCT.

      Materials & methods

      Regularization parameters for the iterative algorithms were found automatically via computer optimization on 60 s acquisitions using the XCAT phantom. Nineteen lung cancer patients were scanned with 60 s arcs using the onboard image on a Varian trilogy linear accelerator. Images were reconstructed using an accelerated ordered subset algorithm. A frequency based temporal regularization algorithm was developed and compared to the McKinnon-Bates algorithm, 4D total variation and prior images compressed sensing (PICCS).

      Results

      All reconstructions were completed in 60 s or less. The proposed method provided a structural similarity of 0.915, compared with 0.786 for the classic McKinnon-bates method. For the patient study, it provided fewer image artefacts than PICCS, and better spatial resolution than 4D TV.

      Conclusion

      Four dimensional iterative CBCT reconstruction was done in less than 60 s, demonstrating the clinical feasibility. The frequency based method outperformed 4D total variation and PICCS on the simulated data, and for patients allowed for tumor location based on 60 s acquisitions, even for slowly breathing patients. It should thus be suitable for routine clinical use.

      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 [
      • 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 [
      • 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 [
      • Wang J.
      • Li T.
      • Xing L.
      Iterative image reconstruction for CBCT using edge-preserving prior.
      ]. For 4D CBCT, it can dramatically reduce scan times [
      • 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 [
      • 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 [
      • Christoffersen C.P.V.
      • Hansen D.
      • Poulsen P.
      • Sørensen T.S.
      Registration-based reconstruction of four-dimensional cone beam computed tomography.
      ,
      • 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.
      ,
      • Zhang Y.
      • Yin F.F.
      • Segars W.P.
      • Ren L.
      A technique for estimating 4D-CBCT using prior knowledge and limited-angle projections.
      ,
      • Mory C.
      • Janssens G.
      • Rit S.
      Motion-aware temporal regularization for improved 4D cone-beam computed tomography.
      ,
      • Liu J.
      • Zhang X.
      • Zhang X.
      • et al.
      5D respiratory motion model based image reconstruction algorithm for 4D cone-beam computed tomography.
      ,
      • Wang J.
      • Gu X.
      High-quality four-dimensional cone-beam CT by deforming prior images.
      ,
      • Wang J.
      • Gu X.
      Simultaneous motion estimation and image reconstruction (SMEIR) for 4D cone-beam CT.
      ,
      • 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 [
      • 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
      minxW(Ax-p)2+R(x)
      (1)


      where A is the projection transform, x is the (4D) image to be reconstructed and p are the log-transformed projections. W is either the identity transform or (as in this work) a diagonal matrix containing weights [
      • Parker D.L.
      Optimal short scan convolution reconstruction for fan beam CT.
      ] to compensate for the offset detectors. Finally, R(x) is the (nonlinear) regularization function, such as total variation, given by
      TV3D(x)=i,j,k,t(xi+1,j,k,t-xi,j,k,t)2+(xi,j+1,k,t-xi,j,k,t)2+(xi,j,k+1,t-xi,j,k,t)2
      (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:
      RF(x)=λTVTV3D(x)+λATVTV3D(Dx)+λF|ReFt(x)|+|ImFt(x)|
      (3)


      D denotes a box filter, downsampling the image by a factor of two in the spatial dimensions. Ft is the Fourier transform along the temporal dimension. λTV, λATV and λF 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. [

      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) [
      • 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 [
      • Christoffersen C.P.V.
      • Hansen D.
      • Poulsen P.
      • Sørensen T.S.
      Registration-based reconstruction of four-dimensional cone beam computed tomography.
      ,
      • 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
      RTV4D(x)=λTVTV3D(x)+λ4DTVt(x)
      (4)


      where TVt(x) is the temporal total variation, given by
      TVt(x)=i,j,k,t|xi,j,k,t+1-xi,j,k,t|
      (5)


      Again, λTV and λ4D are unitless weights, controlling strength of the regularisation.

      2.3 Prior image constrained compressed sensing (PICCS)

      PICCS [
      • 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 xp is reconstructed using the standard Feldkamp-Davis-Kress (FDK) [
      • 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
      RPICCS(x)=λTVTV3D(x)+λPICCStTV3D(xt-xp)
      (6)


      where xt is the 3D image corresponding to phase t of x. The motivation behind PICCS is that most of the image is stationary, so the difference term xt-xp should mainly contain zero. While the original article used a constrained approach Lauzier et al. [
      • 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, λTV and λPICCS 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 [
      • 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 [
      • 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 [

      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 [
      • 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 γtime=0.0002 and γspace=0.0002, 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. [
      • 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 360° 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 [
      • Thirion J.P.
      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 [

      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 [
      • 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 389.12×291.84mm2 with an image resolution of 1024×768pixels. The detector was offset with 144.97 mm, equal to what is used in the clinical scans. 620 projections were obtained over a 360° 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 375×375×375μm3 which was rebinned to 1.5×1.5×1.5mm3.

      2.9 Implementation

      Reconstructions were performed on a NVIDIA Titan X GPU, with 12 GB of RAM. Implementation was done in the CUDA programming language, and is freely available at https://github.com/dchansen/gt-tomography.
      To solve Eq. (1), an accelerated ordered subset (OS) algorithm was used [
      • 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. [
      • 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. [
      • 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. [
      • 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 [
      • Kim D.
      • Ramani S.
      • Fessler J.A.
      Combining ordered subsets and momentum for accelerated X-ray CT image reconstruction.
      ,
      • Nesterov Y.
      Introductory 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 (A and A respectively) is based on our previously published work [
      • 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 [
      • 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 1.5×1.5×1.5mm3. 10 full iterations were chosen, based on previous experience with the same algorithm for other applications [
      • 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 [
      • 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 σ=2.25mm (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 λs’) were found using 500 iterations of surrogate optimization [
      • Mller J.
      • Pich R.
      Mixture surrogate models based on Dempster-Shafer theory for global optimization problems.
      ] using pySOT [

      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.

      3. Results

      3.1 XCAT phantom

      The optimal regularization parameters found for the three iterative methods were of similar magnitude, with PICCS and TV4D having almost the same λTV (Table 1). All the iterative methods had a higher SSIM than FDK and MKB, with the SFR method performs best, closely followed by 4D TV (Table 1). This held true both on average, and for the individual breathing phases (Fig. 1).
      Table 1Reconstruction parameters and resulting SSIM. 4D indicates the regularization in the temporal dimension (λ4D,λPICCS and λF respectively).
      λTVλATVTemporalRuntimeSSIMmin
      TV4D0.20λTV4D=0.6245.3s0.912
      PICCS0.24λPICCS=0.1349.4s0.903
      SFR0.120.42λSFR=0.4956.8s0.916
      FDK0.858
      MKB0.786
      ROOSTER (10 iterations)140s0.857
      ROOSTER1344s0.915
      Figure thumbnail gr1
      Fig. 1SSIM as a function of phase for the different reconstruction methods. Note that the y-axis does not start at 0.
      The ROOSTER method as implemented in RTK required over 20 min computation time to achieve competitive results, but matches the performance of the SFR algorithm. In the reconstructed images the traditional MKB yields a high degree of streak artifacts and the 3D FDK reconstruction exhibits a great deal of blurring, particularly visible at the diaphragm (Fig. 2). Comparing the iterative methods, PICCS has the best contrast in the static soft tissue, where the prior image is reliable. It does however suffer from significant artifacts near the ribs, and the contrast of the tumour is markedly reduced. TV has better tumour contrast, but the relatively thin rib structures are more blurry. SFR retains the good tumour contrast of TV, but has better contrast in the soft tissue region and in the ribs. The cost of this is small artifacts near the diaphragm, which also appear in PICCS and TV, but less markedly. None of the methods are able to fully resolve the blood vessels in the lung.
      Figure thumbnail gr2
      Fig. 2Midventilation phase of the XCAT phantom. Top from left to right: Ground truth, 3D FDK, MKB. Bottom: PICCS, 4DTV and SFR.

      3.2 Patients

      Three representative patients are included in the manuscript. Images from all 19 patients are available in the Supplementary material. Both patient 19 and 15 have a relatively fast breathing period of 2.9 s and 2.3 s respectively, whereas patient 11 has a much slower cycle, of about 6 s (Fig. 3). Patient 11 (Fig. 4) shows the largest difference between the methods. PICCS provides an image with considerable artefacts, and accurate tumour delineation is difficult. TV and SFR are both able to resolve the image clearer. SFR generally displays a sharper image however, and a greater degree of artifact suppression. For patient 15 (Fig. 5), a similar situation can be observed, with PICCS yielding a sharp image with motion artefacts, which cannot be seen with SFR. For patient 19 (shown on Fig. 6), PICCS shows a larger degree of motion artefacts, such as streaking (particularly visible in the axial plane), but has a sharper definition of the static bone structures. TV suffers less from the motion artefacts, but also has a more blurry image. In comparison, SFR demonstrates better resolution, without reintroducing the motion artifacts.
      Figure thumbnail gr3
      Fig. 3Extracted respiration curve for the three patient scans extracted with the method described in Section . Note that the amplitudes extracted do not correlate directly with the diaphragm motion amplitude.
      Figure thumbnail gr4
      Fig. 4Midventilation phase of patient 11 reconstructed using (from left to right) PICCS, 4D TV and SFR.
      Figure thumbnail gr5
      Fig. 5Midventilation phase of patient 15 reconstructed using (from left to right) PICCS, 4D TV and SFR.
      Figure thumbnail gr6
      Fig. 6Midventilation phase of patient 19 reconstructed using (from left to right) PICCS, 4D TV and SFR.
      In comparison with the XCAT phantom, a general trend is that soft tissue contrast is worse in the real scans, where tissue types cannot be distinguished. On the other hand, moving structures like blood vessels are much better resolved.

      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 [
      • 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 [
      • Mory C.
      • Janssens G.
      • Rit S.
      Motion-aware temporal regularization for improved 4D cone-beam computed tomography.
      ,
      • Wang J.
      • Gu X.
      Simultaneous motion estimation and image reconstruction (SMEIR) for 4D cone-beam CT.
      ,
      • Bergner F.
      • Berkus T.
      • Oelhafen M.
      • et al.
      An investigation of 4D cone-beam CT algorithms for slowly rotating scanners.
      ,
      • Vergalasova I.
      • Cai J.
      • Yin F.F.
      A novel technique for markerless, self-sorted 4D-CBCT: feasibility study.
      ,
      • 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. [
      • 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. [
      • 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. [
      • Qi Z.
      • Chen G.H.
      • Rit S.
      • et al.
      Performance studies of four-dimensional cone beam computed tomography.
      ] and Gao et al. [
      • 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. [
      • 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. [
      • 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.

      Conflict of interest

      We have no conflicts of interest to declare.

      Acknowledgements

      The authors would like to thank Mai Lykkegaard Schmidt for her help in collecting the clinical data. This work was supported by the Danish Cancer Society (Grant No. R90-A5992). Thomas Sangild Sørensen is also employed at 3Shape Inc.

      Appendix – Algorithm

      Tabled 1
      Algorithm 1. Momentum accelerated ordered subsets method with augmented Lagrangian multipliers for penalized weighted least squares (PWLS).
      1
      2 function OS-MOM(f)
      3  x00
      4  gAWWA1    ▷SQS diagonal preconditioner
      5  t01
      6  for k1 to iterations do
      7   for l1 to subsets do   ▷Loop over subsets
      8    ml+k×subsets
      9    ymsubsets·Al(fl-Alxm-1)øgø: elementwise-division
      10    ym denoise(ym,g)   ▷Denoise using algorithm 2
      11   tm121+1+4tm-1·tm-1
      12   xmym+tm-1-1tmym-ym-1 ▷ Momentum acceleration step
      13  return xfinal
      Tabled 1
      Algorithm 2. Arrow Hurwicz denoising
      • Arrow K.J.
      • Hurwicz L.
      • Uzawa H.
      Studies in linear and non-linear programming.
      ,
      • Chambolle A.
      • Pock T.
      A first-order primal-dual algorithm for convex problems with applications to imaging.
      .
      1 function denoise(x0,g)
      2  τ110-3
      3  γ0.35
      4  σ1116τ1
      5  for k1 to iterations do
      6   pkyk-1+σkKx¯k-1   ▷ K calculates the finite difference
      7   yikpik1+σkαmax1,pik1+σkα
      8   xkxk-1-τkKyk+τkλx0øg
      9   xkxkøg
      10   θk11+2γτ
      11   τk+1θkτk
      12   σk+1σkθk
      13  return P(xfinal)   ▷ P sets all negative elements to 0

      Supplementary data

      References

        • 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.
        Radiat Oncol. 2012; 7: 81
        • 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.
        J Appl Clin Med Phys. 2016; 17: 97-106
        • Wang J.
        • Li T.
        • Xing L.
        Iterative image reconstruction for CBCT using edge-preserving prior.
        Med Phys. 2008; 36: 252-260
        • Qi Z.
        • Chen G.H.
        • Rit S.
        • et al.
        Performance studies of four-dimensional cone beam computed tomography.
        Phys Med Biol. 2011; 56: 6709-6721
        • Christoffersen C.P.V.
        • Hansen D.
        • Poulsen P.
        • Sørensen T.S.
        Registration-based reconstruction of four-dimensional cone beam computed tomography.
        IEEE Trans Med Imaging. 2013; 32: 2064-2077
        • 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.
        Med Phys. 2017; 44: 1089-1104
        • Zhang Y.
        • Yin F.F.
        • Segars W.P.
        • Ren L.
        A technique for estimating 4D-CBCT using prior knowledge and limited-angle projections.
        Med Phys. 2013; 40: 121701
        • Mory C.
        • Janssens G.
        • Rit S.
        Motion-aware temporal regularization for improved 4D cone-beam computed tomography.
        Phys Med Biol. 2016; 61: 6856-6877
        • Liu J.
        • Zhang X.
        • Zhang X.
        • et al.
        5D respiratory motion model based image reconstruction algorithm for 4D cone-beam computed tomography.
        Inverse Prob. 2015; 31: 115007
        • Wang J.
        • Gu X.
        High-quality four-dimensional cone-beam CT by deforming prior images.
        Phys Med Biol. 2013; 58: 231-246
        • Wang J.
        • Gu X.
        Simultaneous motion estimation and image reconstruction (SMEIR) for 4D cone-beam CT.
        Med Phys. 2013; 40: 101912
        • Zhong Z.
        • Gu X.
        • Mao W.
        • Wang J.
        4D cone-beam CT reconstruction using multi-organ meshes for sliding motion modeling.
        Phys Med Biol. 2016; 61: 996-1020
        • Parker D.L.
        Optimal short scan convolution reconstruction for fan beam CT.
        Med Phys. 1982; 9: 254-257
      1. 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.

        • Gao H.
        • Li R.
        • Lin Y.
        • Xing L.
        4D cone beam CT via spatiotemporal tensor framelet.
        Med Phys. 2012; 39: 6943-6946
        • 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.
        Med Phys. 2008; 35: 660-663
        • Feldkamp L.A.
        • Davis L.C.
        • Kress J.W.
        Practical cone-beam algorithm.
        J Opt Soc Am A. 1984; 1: 612
        • Lauzier P.T.
        • Tang J.
        • Chen G.H.
        Prior image constrained compressed sensing: implementation and performance evaluation.
        Med Phys. 2011; 39: 66-80
        • Mc Kinnon G.C.
        • Bates R.H.
        Towards imaging the beating heart usefully with a conventional CT scanner.
        IEEE Trans Bio-med Eng. 1981; 28: 123-127
        • 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).
        J Phys Conf Ser. 2014; 489 (012079)
      2. 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.

        • 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.
        Acta Oncol. 2014; 53: 1107-1113
        • Thirion J.P.
        Image matching as a diffusion process: an analogy with Maxwell’s demons.
        Med Image Anal. 1998; 2: 243-260
      3. 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.

        • Segars W.P.
        • Sturgeon G.
        • Mendonca S.
        • Grimes J.
        • Tsui B.M.W.
        4D XCAT phantom for multimodality imaging research.
        Med Phys. 2010; 37: 4902
        • Hansen D.C.
        • Sangild Sørensen T.
        • Rit S.
        Fast reconstruction of low dose proton CT by sinogram interpolation.
        Phys Med Biol. 2016; 61: 5868-5882
        • Kim D.
        • Ramani S.
        • Fessler J.A.
        Combining ordered subsets and momentum for accelerated X-ray CT image reconstruction.
        IEEE Trans Med Imaging. 2015; 34: 167-178
        • Arrow K.J.
        • Hurwicz L.
        • Uzawa H.
        Studies in linear and non-linear programming.
        in: Stanford Mathematical Studies in the Social Sciences. vol. II. Stanford University Press, Stanford, California1958
        • 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.
        Med Phys. 2016; 43: 1849-1872
        • Nesterov Y.
        Introductory lectures on convex optimization: a basic course.
        Springer US, New York2004
        • Wang Z.
        • Bovik A.
        • Sheikh H.
        • Simoncelli E.
        Image quality assessment: from error visibility to structural similarity.
        IEEE Trans Image Process. 2004; 13: 600-612
        • Mller J.
        • Pich R.
        Mixture surrogate models based on Dempster-Shafer theory for global optimization problems.
        J Global Optim. 2011; 51: 79-104
      4. Eriksson D, Bindel D, Shoemaker C. Surrogate Optimization Toolbox (pySOT) [Internet]. Cited [2018 Mar 2]. Available from: github.com/dme65/pySOT.

        • Bergner F.
        • Berkus T.
        • Oelhafen M.
        • et al.
        An investigation of 4D cone-beam CT algorithms for slowly rotating scanners.
        Med Phys. 2010; 37: 5044
        • Vergalasova I.
        • Cai J.
        • Yin F.F.
        A novel technique for markerless, self-sorted 4D-CBCT: feasibility study.
        Med Phys. 2012; 39: 1442
        • 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.
        Med Phys. 2014; 41 (041912)
        • Yan H.
        • Zhen X.
        • Folkerts M.
        • et al.
        A hybrid reconstruction algorithm for fast and accurate 4D cone-beam CT imaging.
        Med Phys. 2014; 41 (071903)
        • Chambolle A.
        • Pock T.
        A first-order primal-dual algorithm for convex problems with applications to imaging.
        J Math Imaging Vision. 2011; 40: 120-145