Fast Statistical Iterative Reconstruction for MVCT in TomoTherapy

TL;DR

This paper presents a fast statistical iterative reconstruction algorithm for MVCT in TomoTherapy, significantly improving image quality and speed using GPU parallelization and optimized iterative techniques.

Contribution

The authors developed a GPU-accelerated iterative reconstruction method that produces high-quality MVCT images within seconds, outperforming traditional FBP in both speed and image quality.

Findings

SNR and CNR were several times better than FBP.

GPU parallelization achieved several hundred times faster reconstruction.

High-quality images obtained within a few seconds.

Abstract

Statistical iterative reconstruction is expected to improve the image quality of megavoltage computed tomography (MVCT). However, one of the challenges of iterative reconstruction is its large computational cost. The purpose of this work is to develop a fast iterative reconstruction algorithm by combining several iterative techniques and by optimizing reconstruction parameters. Megavolt projection data was acquired from a TomoTherapy system and reconstructed using our statistical iterative reconstruction. Total variation was used as the regularization term and the weight of the regularization term was determined by evaluating signal-to-noise ratio (SNR), contrast-to-noise ratio (CNR), and visual assessment of spatial resolution using Gammex and Cheese phantoms. Gradient decent with an adaptive convergence parameter, ordered subset expectation maximization (OSEM), and CPU/GPU parallelization were applied in order to accelerate the present reconstruction algorithm. The SNR and CNR of the iterative reconstruction were several times better than that of filtered back projection (FBP). The GPU parallelization code combined with the OSEM algorithm reconstructed an image several hundred times faster than a CPU calculation. With 500 iterations, which provided good convergence, our method produced a 512$\times$512 pixel image within a few seconds. The image quality of the present algorithm was much better than that of FBP for patient data. An image from the iterative reconstruction in TomoTherapy can be obtained within few seconds by fine-tuning the parameters. The iterative reconstruction with GPU was fast enough for clinical use, and largely improve the MVCT images.

Figures

Tables

Table 1. Table 1: Reconstruction speed [s] per one slice image with 512 × 512 512 512 512\times 512 .

Thinning parameter	1	1/2	1/5	1/10	1/20	1/40
CPU (12 threads)	442.91	244.34	138.77	80.68	58.79	43.29
CPU (20 threads)	325.40	191.16	94.40	51.21	35.99	24.59
GPU (GTX1080Ti)	11.00	7.76	5.37	4.30	3.71	3.39

Equations

$${\mathbf{\mu }}^{\ast }$$

$$\ln P({\mathbf{\mu }}^{\ast }|\mathbf{y})\sim \ln P(\mathbf{y}|{\mathbf{\mu }}^{\ast })+\ln P({\mathbf{\mu }}^{\ast })=L({\mathbf{\mu }}^{\ast })+\alpha R({\mathbf{\mu }}^{\ast }).$$

$${\mu }_r^{\ast (n+1)}={\mu }_r^{\ast (n)}+{\lambda }^{(n)}\frac{\partial }{\partial {\mu }_r^{\ast (n)}}\left\{L({\mathbf{\mu }}^{\ast (n)})+\alpha R({\mathbf{\mu }}^{\ast (n)})\right\},$$

$${\lambda }^{(n)}={\lambda }^{(n-1)}+\mathrm{0.0001.}$$

$${\lambda }^{(n)}={\lambda }^{(n-1)}-\mathrm{0.0002.}$$

Full text

Technical Note: Fast Statistical Iterative Reconstruction for MVCT in TomoTherapy

Sho Ozaki1, Akihiro Haga2a, Edward Chao3, Calvin Maurer3, Kanabu Nawa1, Takeshi Ohta1, Takahiro Nakamoto1, Yuki Nozawa1, Taiki Magome4, Masahiro Nakano5, Keiichi Nakagawa1

1*Department of Radiology, University of Tokyo Hospital, Tokyo, Japan,

2Graduate School of Biomedical Sciences, Tokushima University, Tokushima, Japan*

3*Accuray Incorporated. Sunnyvale, CA, USA

4Radiological Science, Komazawa University, Tokyo, Japan

5Radiation Oncology Department, Cancer Institute Hospital, Japanese Foundation for Cancer Research, Tokyo, Japan*

Numbers of Text pages, Figures, and Tables

Total pages: 11 (including title page), Tables: 1, Figures: 4

Running head

Fast Iterative Reconstruction in TomoTherapy

Contact information

a) Electronic mail: haga@tokushima-u.ac.jp

Telephone number: +81 88 633 9024

Fax number: +81 88 633 9024

Mailing address: 3-18-15, Kuramoto, Tokushima, 770-8503, Japan

Abstract

Purpose: Statistical iterative reconstruction is expected to improve the image quality of megavoltage computed tomography (MVCT). However, one of the challenges of iterative reconstruction is its large computational cost. The purpose of this work is to develop a fast iterative reconstruction algorithm by combining several iterative techniques and by optimizing reconstruction parameters.

Methods: Megavolt projection data was acquired from a TomoTherapy system and reconstructed using our statistical iterative reconstruction. Total variation was used as the regularization term and the weight of the regularization term was determined by evaluating signal-to-noise ratio (SNR), contrast-to-noise ratio (CNR), and visual assessment of spatial resolution using Gammex and Cheese phantoms. Gradient decent with an adaptive convergence parameter, ordered subset expectation maximization (OSEM), and CPU/GPU parallelization were applied in order to accelerate the present reconstruction algorithm.

Results: The SNR and CNR of the iterative reconstruction were several times better than that of filtered back projection (FBP). The GPU parallelization code combined with the OSEM algorithm reconstructed an image several hundred times faster than a CPU calculation. With 500 iterations, which provided good convergence, our method produced a 512$\times$512 pixel image within a few seconds. The image quality of the present algorithm was much better than that of FBP for patient data.

Conclusion: An image from the iterative reconstruction in TomoTherapy can be obtained within few seconds by fine-tuning the parameters. The iterative reconstruction with GPU was fast enough for clinical use, and largely improve the MVCT images.

I Introduction

Helical TomoTherapy (HT) (Accuray, Sunnyvale, CA) is an innovative machine that delivers intensity modulated radiation therapy (IMRT) Mackie:1993 ; Mackie:2006 , and includes megavoltage computed tomography (MVCT) for image guided radiation therapy (IGRT) Ruchala:1999 ; Meeks:2005 . In order to ensure precise image guidance, the image quality of the MVCT is crucial. In particular, detection of soft tissue contrast is very important for soft tissue based registration accuracy and the accuracy of adaptive radiotherapy depends on the accuracy of the soft tissue deformation during the daily treatment Langen:2005 . However, in general, MVCT images are noisier than that of kilovoltage CT (kVCT) Ruchala:1999 . Furthermore, the relative difference in attenuation coefficients of various soft tissues at MV energies is smaller than at kV energies. The smaller soft tissue attenuation differences lead to smaller contrast differences in MVCT and make it more difficult to register soft tissues accurately.

Statistical iterative reconstruction can improve the image quality of MVCT without increasing the imaging dose. Unlike conventional reconstruction methods such as the filtered back projection (FBP), the iterative reconstruction methods can take into account a priori information including the CT geometry, photon statistics at the detector, and known properties of the image, in addition to the observed projection data. With such knowledge-based information, the iterative reconstruction method can reduce noise and improve soft tissue contrast. Accuray has introduced an iterative reconstruction algorithm, ”CTrue IR”, on their Radixact platform that appears to improve the MVCT image quality AccurayURL . However, there aren’t any publications describing the Accuray’s algorithm as well as any source codes available. This is one of the motivations of this paper despite the availability of the Accuray CTrue IR.

A disadvantage of the iterative reconstruction is the large computational cost. Since the method iteratively calculates a reprojection to compare with measured projection data, the total calculation time can be large. In this note, we describe an acceleration of the statistical iterative reconstruction for TomoTherapy MVCT by utilizing a parallelized CPU/GPU calculation. There are several reports on GPU based iterative reconstruction for cone-beam CT (CBCT)  Scherl: 2007 ; Jia:2010 ; Xu:2012 . The reconstruction time for a set of CBCT images with the iterative reconstruction algorithm could be shortened to $77$ to $130$s (about $100$ times faster than with conventional approaches). In addition to the GPU parallelization, we applied several acceleration techniques including an adaptive convergence parameter for the gradient decent method and the ordered subset expectation maximization (OSEM) algorithm. The performance and the quality of the reconstructed images of our method were evaluated using phantom and patient data.

II Iterative reconstruction scheme

In this section, we briefly explain the basics of the statistical iterative CT reconstruction. In this study, the maximum a posteriori (MAP) approach was employed as our image reconstruction frameworkTang2009a ; Hanson1983a ; Green1990 ; Lange1995 . The concept of the MAP approach is to maximize the posteriori probability of reconstruction images with observed projection (or sinogram) data and a prior information about the image, encoded in the regularization term. Reconstructed images can be obtained from an iterative process to maximize the log a posteriori probability function $\ln P(\bm{\mu}^{*}|\bm{y})$ as,

[TABLE]

where

[TABLE]

Here $L({\bm{\mu}}^{*})=\ln P(\bm{y}|\bm{\mu}^{*})$ expresses the log-likelihood probability function of observing the projection data set $\bm{y}$ at the given expectation of the image $\bm{\mu}^{*}$, whereas $\alpha R({\bm{\mu}}^{*})=\ln P(\bm{\mu}^{*})$ stands for the regularization term in the optimization process with a hyper parameter $\alpha$ controlling the weight of the regularization. We assumed that the photon statistics at the detector obeys the Poisson distribution and employed the total variation (TV) as the regularization term. The TV term penalizes a large difference between the image values of a certain position and its neighbors, and thus reduces the noise of the image. At the same time, however, this term also reduces the sharpness of object edges in the image, which leads to edge blurring. Therefore, in general, the benefit of noise reduction must be balanced against the negative impact of edge blurring. We seek the appropriate value of the weight $\alpha$ by examining the trade-off.

In the optimization of Eq. (1), the gradient decent method including the parameter adaptation was employed,

[TABLE]

where $n$ is the number of the iteration step, and $r$ means the index of the pixel to be reconstructed. $\lambda^{(n)}$ is the parameter which is adjusted to accelerate the calculation iteration step by iteration step as described in the next section. It is noted that, as well as the value of $\lambda^{(n)}$, the initial image setting $\mu_{r}^{*(0)}$ can affect the convergence speed.

III Setups and rapid convergence methods

III.1 Data acquisition and reconstruction

Gammex phantom (Gammex, Middleton, WI) and Cheese phantom (Accuray, Sunnyvale, CA) were used in the analysis as well as the two patients for visual demonstration. All the data were acquired at the University of Tokyo Hospital during 1st - 31st May, 2016. The reconstruction size of the image is $512\times 512$ with an equal pixel size of 1mm for the one slice.

III.2 Computer used in reconstruction

The CPU/GPU specs in the computer used in this study (HP customized workstation: Z840/Z800) are shown as follows:

• •

CPU: Xeon(R) X5690 ($\times 2$: $12$ threads) / Xeon(R) E5-2687Wv ($\times 2$: $20$ threads)

• •

GPU: GeForce GTX 1080 Ti

III.3 Parameter optimization

III.3.1 Total variation parameter

The weight of the TV regularization term $\alpha$ was set as [math], $0.0001$, $0.0002$, and $0.0003$. The optimal values of the weight parameter were determined by estimating the signal-to-noise ratio (SNR), the contrast-to-noise ratio (CNR), the visual resolution, and the edge blurring effect. For SNR/CNR and edge blurring analyses, the Gammex phantom was used, whereas for visual resolution, the Cheese phantom was used. The region-of-interests (ROIs) used in the SNR/CNR analysis as well as their definitions and the edge blurring analysis method are mentioned in the supplemental material.

III.3.2 Convergence parameter and initial image setting

Although a relatively large value of $\lambda^{(n)}$ in Eq.(3) gives a rapid convergence, a too large value leads to non-convergence. Also the best value could depend on the iteration step. Therefore, we employed the following adaptive convergence parameter which depends on the iteration step: when the objective function being the square of difference between a reprojection and a projection decreases, the convergence parameter is increased as

[TABLE]

When the objective function increases, the parameter is decreased as

[TABLE]

The initial value of this parameter was set as 0.01.

As mentioned, the initial image setting $\mu_{r}^{*(0)}$ can also affect the convergence speed. In this study, homogeneous images with $\mu_{r}^{*(0)}=0$, $0.05$, and $0.1$ [cm*-1*] were tested.

III.4 Ordered subset expectation maximization

The ordered subset expectation maximization (OSEM) is a well-known method to accelerate the convergence of the iterative reconstruction algorithm Hudson:1994 . A basic idea of OSEM is to divide the original projection data into a subset of the data and to use the different subset in each iteration step. The easiest way to divide the projection data into subset is to thin the sinogram height out. In the HT system, the data sampling time is $80$ Hz, and it takes $10$ sec for one rotation. Hence the sinogram height is $800$. We thinned the sinogram height out by $1/2$, $1/5$, $1/10$, $1/20$, and $1/40$ (defined as the thinning parameter) with an equal interval (the sinogram height used in each iteration step is $400$, $160$, $80$, $40$, and $20$, respectively). For instance, in the $1/40$ case, $40$ iterations are needed to cover the whole sinogram data, so that the convergence could be relatively slow. Therefore, the convergence behavior as well as the convergence efficiency has to be verified carefully. Here, it should be noted that the regularization weight $\alpha$ in the OSEM method is multiplied by the thinning parameter in order to keep the ratio of log-likelihood function with regularization term.

III.5 Parallel computing

For the parallel computing, we applied openMP for CPU parallelization, and cuda for GPU parallelization. We have made our source code available at https:$\cdots$. In the next section, the reconstruction speeds with the several cases of the number of thinning out for OSEM described above will be compared.

IV Results and discussion

IV.1 Total variation parameter

In order to determine the total variation parameter, the SNR and the CNR in the eleven insertions of the Gammex phantom were analyzed (see Fig. S-2 in the supplemental material for the reconstructed images). Figure 1 shows the results of the SNR and the CNR. In these results, the SNR and the CNR of the iterative reconstruction are normalized by the values of the FBP. In Fig. 1, both of the SNR and the CNR converge as the iteration number increases, and their values with the finite $\alpha$ are always larger than that of the FBP image. On the other hand, as mentioned above, noise suppression is generally accompanied by a trade-off of reduced image sharpness. This is seen in Fig. 2, where we can see that as the total variation parameter increases, the noise of the image decreases but the visibility of the small holes in the resolution plug of the cheese phantom also decreases. From these results, we concluded that $0.0001\leq\alpha\leq 0.0002$ maintains acceptable image resolution while achieving significant noise suppression.

IV.2 Convergence parameter and initial image setting

Figure 3 shows the convergence behavior for a variety of parameters, where convergence is characterized by the square of the difference between a projection and a reprojection in the sinogram domain. As expected, better convergence is obtained with the adaptive convergence parameter method. With the adaptive convergence method, a similar degree of convergence is achieved with 500 iterations compared to 1500 iterations using a constant convergence parameter. It should be noted that the oscillatory behavior of the square difference is also caused by this adaptation; As the convergence parameter becomes gradually large according to the Eq. (4), the images grow away from the optimal values. Eventually, the square difference turns to increase and numerically diverge. This increasing behavior is however suppressed by decreasing the parameter as Eq. (5) when the square difference is increased. Due to this oscillation, the reconstructed image at the minimum of the oscillatory behavior of the square difference was employed. The initial homogeneous images were set to values of [math], $0.05$ and $0.1$ [cm*-1*], and the best convergence was obtained for $0.05$ [cm*-1*] in the Gammex phantom, though the difference was so small. This was an expected result because the reconstructed value of $\mu^{*}$ in the water area was approximately 0.07 [cm*-1*] and the mean value of the reconstructed image including the air area was 0.026 [cm*-1*]. Other two phantoms with the different radius also gave the same result (see Table S-1 in supplemental file).

It was found that the adaptive convergence parameter and the initial image setting accelerated the convergence approximately three times faster than that with the constant convergence parameter starting with the “zero” image.

IV.3 CPU parallelization and OSEM method

The convergence behavior with applying the OSEM method is also shown in Fig. 3, where the oscillation behavior is emphasized. The convergence needs relatively larger iteration steps in the OSEM method, because the reconstructed image and its reprojection are obtained from thinned-out projection data. However, this apparent disadvantage is overcome by large reduction of the calculation time for the one iteration, thanks to the thinning. Eventually, the OSEM method highly accelerates the reconstruction speed. We found that 500 iterations were an enough number of iteration in which the square of difference converges in the case of the thinning parameter more than $1/20$, whereas more iterations was required in the thinning parameter of $1/40$. Combining OSEM with CPU parallelization yielded the image reconstruction time more than $40$ times faster than the original single CPU implementation, and more than $8$ times faster than the CPU parallelization without OSEM. Consequently, combining the OSEM with CPU parallelization considerably accelerates the calculation.

IV.4 GPU computing

In our reconstruction algorithm, the reprojection and the estimation processes occupy a large fraction of the calculation time. Therefore, we coded these parts of the reconstruction using CUDA GPU processing. Table I shows the average calculation time for reconstructing one image slice. Although the results depend on the iteration number and OSEM parameter, the GPU code produces an image more than $10$ times faster than the parallelized CPU code without degrading image quality. Using the GPU code, a high quality image can be reconstructed at a speed of $3.71$ seconds per slice. This large improvement in iterative reconstruction speed on the GPU is achieved without any degradation in image quality compared to the CPU reconstructed images.

IV.5 Reconstructed images

Finally, the representative reconstructed images with $\alpha=0.0002$ are shown in Fig. 4. For comparison, FBP reconstruction images are also shown there. In addition to the noise reduction seen in Figs. 4(a)-4(b), a relatively large reconstructed region (field of view) was achieved with the iterative reconstruction approach seen in Figs 4(c)-(d). Owing to the fan beam geometry, information of the projection data decrease at a fringe region of FOV. Therefore, FBP fails to reconstruct the image at the edge of the body (Fig. 4-(c)). On the other hand, the statistical iterative reconstruction algorithm can reconstruct the image with a fewer projection data (the examination of the reconstruction using a half of projection data is given in supplemental file), and thus shows much better image quality at the edge (Fig. 4-(d). From the viewpoint of the image quality, it can be said that there is a large benefit to use the iterative reconstruction algorithm clinically.

V Conclusion

We investigated the clinical feasibility of a statistical iterative reconstruction method for MVCT on a HT system. An image from the iterative reconstruction can be obtained within few seconds per image slice with a size of $512\times 512$. In addition, the image quality of our iterative method is much better than that of the FBP. It can be concluded that the iterative reconstruction with GPU is fast enough for clinical use, and largely improves the MVCT images.

**Conflict of interest

**Edward Chao and Calvin Maurer are officers of Accuray Inc.

**Funding

**Keiichi Nakagawa has a research grant from Accuray Inc.

**Ethical approval

**The study was ethically approved by the institutional review board at the University of Tokyo Hospital (reference number 3372).