Fast 4D cone-beam CT from 60 s acquisitions.

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.

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]. 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]. Iterative cone-beam CT reconstruction has been shown to allow for reduced imaging dose and improved image quality [3]. For 4D CBCT, it can dramatically reduce scan times [4]. On the other hand, reconstruction times are quite long ranging from 5 min to several hours in literature [5]. 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], [6], [7], [8], [9], [10], [11], [12]. 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]. 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.

Materials and methods

Iterative CBCT reconstruction is often treated as a regularized least squares problemwhere 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] to compensate for the offset detectors. Finally, is the (nonlinear) regularization function, such as total variation, given byin 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.

Sparse Frequency Regularization (SFR)

The following regulariser is proposed: 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] 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], but only the first-order differential was included, in the interest of retaining clinical reconstruction times.

4D total variation

Several previous publications have demonstrated the effective use of 4D TV for CT reconstruction [5], [8]. 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 bywhere is the temporal total variation, given byAgain, and are unitless weights, controlling strength of the regularisation.

Prior image constrained compressed sensing (PICCS)

PICCS [16] is an extension of standard 3D TV. First, a 3D image is reconstructed using the standard Feldkamp-Davis-Kress (FDK) [17] method. This is then used as a prior, and the regularization term becomeswhere 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] 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.

McKinnon-Bates

The McKinnon-Bates (MKB) algorithm [19] 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.

RTK ROOSTER

The RTK framework [20] 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] implemented in this framework. We do not use the motion-aware version [8], 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.

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]. 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.

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]. 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], as the diaphragm was not consistently visible in all projections causing that method to fail.

Virtual XCAT phantom

To validate our approach and tune the regularisation parameters the female XCAT 4D phantom [25] 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 .

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 [26] (see Appendix). This is a modification of the algorithm proposed by Kim et al. [27], adapted to allow for isotropic total variation by combining it with the work of Arrow et al. [28]. The approach is similar to the one recently presented by Xu et al. [29] 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.

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], [30] 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.

Forward and backward projection

The GPU implementation of the basic forward and backward projection ( and respectively) is based on our previously published work [5], 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].

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]. To evaluate image quality on the XCAT phantom, structural similarity (SSIM) was used [31]. 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] using pySOT [33]. The found parameters were then used for the patient scans.

Results

XCAT phantom

The optimal regularization parameters found for the three iterative methods were of similar magnitude, with PICCS and TV4D having almost the same (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 1

Reconstruction parameters and resulting SSIM. 4D indicates the regularization in the temporal dimension ( and respectively).

	λ_TV	λ_ATV	Temporal	Runtime	SSIM_min
TV4D	0.20	–	λ_TV4D=0.62	45.3s	0.912
PICCS	0.24	–	λ_PICCS=0.13	49.4s	0.903
SFR	0.12	0.42	λ_SFR=0.49	56.8s	0.916
FDK	–	–	–	–	0.858
MKB	–	–	–	–	0.786
ROOSTER (10 iterations)	–	–	–	140s	0.857
ROOSTER	–	–	–	1344s	0.915

Fig. 1

SSIM 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.

Fig. 2

Midventilation phase of the XCAT phantom. Top from left to right: Ground truth, 3D FDK, MKB. Bottom: PICCS, 4DTV and SFR.

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.

Fig. 3

Extracted respiration curve for the three patient scans extracted with the method described in Section 2.7. Note that the amplitudes extracted do not correlate directly with the diaphragm motion amplitude.

Fig. 4

Midventilation phase of patient 11 reconstructed using (from left to right) PICCS, 4D TV and SFR.

Fig. 5

Midventilation phase of patient 15 reconstructed using (from left to right) PICCS, 4D TV and SFR.

Fig. 6

Midventilation 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.

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], 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], [11], [34], [35], [36].

A few other works have investigated 4D CBCT from 60 s acquisitions such as Yan et al. [37] and Zhong et al. [12]. 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] and Gao et al. [15] 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] investigated a new registration based algorithm and Schmidt et al. [22] 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.

Algorithm 1. Momentum accelerated ordered subsets method with augmented Lagrangian multipliers for penalized weighted least squares (PWLS).
1
2 function OS-MOM(f)
3  x^0←0
4  g←A^†W^†WA1    ▷SQS diagonal preconditioner
5  t^0←1
6  fork←1 to iterations do
7   forl←1 to subsets do   ▷Loop over subsets
8    m←l+k×subsets
9    y^m∗←subsets·Al†(f_l-A_lx^m-1)øg ▷ ø: elementwise-division
10    y^mdenoise(y^m∗,g)   ▷Denoise using algorithm 2
11   t^m←121+1+4t^m-1·t^m-1
12   x^m←y^m+t^m-1-1t^my^m-y^m-1 ▷ Momentum acceleration step
13  returnx^final

Algorithm 2. Arrow Hurwicz denoising [28], [38].
1 functiondenoise(x^0,g)
2  τ^1←10^-3
3  γ←0.35
4  σ^1←116τ^1
5  fork←1 to iterations do
6   p^k←y^k-1+σ^kKx¯^k-1   ▷ K calculates the finite difference
7   yik←pik1+σ^kαmax1,pik1+σ^kα
8   x^k∗←x^k-1-τ^kK^†y^k+τ^kλx^0øg
9   x^k←x^k∗øg
10   θ^k←11+2γτ
11   τ^k+1←θ^kτ^k
12   σ^k+1←σ^kθ^k
13  returnP(x^final)   ▷ P sets all negative elements to 0