Search

This article examines the critical role of fast Monte Carlo dose calculations in advancing proton therapy techniques, particularly in the context of increasing treatment customization and precision. As adaptive radiotherapy and other patient-specific approaches evolve, the need for accurate and precise dose calculations, essential for techniques like proton-based stereotactic radiosurgery, becomes more prominent. These calculations, however, are time-intensive, with the treatment planning/optimization process constrained by the achievable speed of dose computations. Thus, enhancing the speed of Monte Carlo methods is vital, as it not only facilitates the implementation of novel treatment modalities but also improves the optimality of treatment plans. Today, the state-of-the-art in Monte Carlo dose calculation speeds is 106 - 107 protons per second. This review highlights the latest advancements in fast Monte Carlo dose calculations that have led to such speeds, including emerging artificial intelligence-based techniques, and discusses their application in both current and emerging proton therapy strategies., Comment: 41 pages

Inverse treatment planning in radiation therapy is formulated as optimization problems. The objective function and constraints consist of multiple terms designed for different clinical and practical considerations. Weighting factors of these terms are needed to define the optimization problem. While a treatment planning system can solve the optimization problem with given weights, adjusting the weights for high plan quality is performed by human. The weight tuning task is labor intensive, time consuming, and it critically affects the final plan quality. An automatic weight-tuning approach is strongly desired. The weight tuning procedure is essentially a decision making problem. Motivated by the tremendous success in deep learning for decision making with human-level intelligence, we propose a novel framework to tune the weights in a human-like manner. Using treatment planning in high-dose-rate brachytherapy as an example, we develop a weight tuning policy network (WTPN) that observes dose volume histograms of a plan and outputs an action to adjust organ weights, similar to the behaviors of a human planner. We train the WTPN via end-to-end deep reinforcement learning. Experience replay is performed with the epsilon greedy algorithm. Then we apply the trained WTPN to guide treatment planning of testing patient cases. The trained WTPN successfully learns the treatment planning goals to guide the weight tuning process. On average, the quality score of plans generated under the WTPN's guidance is improved by ~8.5% compared to the initial plan with arbitrary weights, and by 10.7% compared to the plans generated by human planners. To our knowledge, this is the first tool to adjust weights for the treatment planning in a human-like fashion based on learnt intelligence. The study demonstrates potential feasibility to develop intelligent treatment planning system via deep reinforcement learning.

A number of image-processing problems can be formulated as optimization problems. The objective function typically contains several terms specifically designed for different purposes. Parameters in front of these terms are used to control the relative weights among them. It is of critical importance to tune these parameters, as quality of the solution depends on their values. Tuning parameter is a relatively straightforward task for a human, as one can intelligently determine the direction of parameter adjustment based on the solution quality. Yet manual parameter tuning is not only tedious in many cases, but becomes impractical when a number of parameters exist in a problem. Aiming at solving this problem, this paper proposes an approach that employs deep reinforcement learning to train a system that can automatically adjust parameters in a human-like manner. We demonstrate our idea in an example problem of optimization-based iterative CT reconstruction with a pixel-wise total-variation regularization term. We set up a parameter tuning policy network (PTPN), which maps an CT image patch to an output that specifies the direction and amplitude by which the parameter at the patch center is adjusted. We train the PTPN via an end-to-end reinforcement learning procedure. We demonstrate that under the guidance of the trained PTPN for parameter tuning at each pixel, reconstructed CT images attain quality similar or better than in those reconstructed with manually tuned parameters., Comment: 8 pages, 8 figures, 2 tables

We recently built an analytical source model for GPU-based MC dose engine. In this paper, we present a sampling strategy to efficiently utilize this source model in GPU-based dose calculation. Our source model was based on a concept of phase-space-ring (PSR). This ring structure makes it effective to account for beam rotational symmetry, but not suitable for dose calculations due to rectangular jaw settings. Hence, we first convert PSR source model to its phase-space let (PSL) representation. Then in dose calculation, different types of sub-sources were separately sampled. Source sampling and particle transport were iterated. So that the particles being sampled and transported simultaneously are of same type and close in energy to alleviate GPU thread divergence. We also present an automatic commissioning approach to adjust the model for a good representation of a clinical linear accelerator . Weighting factors were introduced to adjust relative weights of PSRs, determined by solving a quadratic minimization problem with a non-negativity constraint. We tested the efficiency gain of our model over a previous source model using PSL files. The efficiency was improved by 1.70 ~ 4.41, due to the avoidance of long data reading and transferring. The commissioning problem can be solved in ~20 sec. Its efficacy was tested by comparing the doses computed using the commissioned model and the uncommissioned one, with measurements in different open fields in a water phantom under a clinical Varian Truebeam 6MV beam. For the depth dose curves, the average distance-to-agreement was improved from 0.04~0.28 cm to 0.04~0.12 cm for build-up region and the root-mean-square (RMS) dose difference after build-up region was reduced from 0.32%~0.67% to 0.21%~0.48%. For lateral dose profiles, RMS difference was reduced from 0.31%~2.0% to 0.06%~0.78% at inner beam and from 0.20%~1.25% to 0.10%~0.51% at outer beam., Comment: page line has been removed

VMAT optimization is a computationally challenging problem due to its large data size, high degrees of freedom, and many hardware constraints. High-performance graphics processing units have been used to speed up the computations. However, its small memory size cannot handle cases with a large dose-deposition coefficient (DDC) matrix. This paper is to report an implementation of our column-generation based VMAT algorithm on a multi-GPU platform to solve the memory limitation problem. The column-generation approach generates apertures sequentially by solving a pricing problem (PP) and a master problem (MP) iteratively. The DDC matrix is split into four sub-matrices according to beam angles, stored on four GPUs in compressed sparse row format. Computation of beamlet price is accomplished using multi-GPU. While the remaining steps of PP and MP problems are implemented on a single GPU due to their modest computational loads. A H&N patient case was used to validate our method. We compare our multi-GPU implementation with three single GPU implementation strategies: truncating DDC matrix (S1), repeatedly transferring DDC matrix between CPU and GPU (S2), and porting computations involving DDC matrix to CPU (S3). Two more H&N patient cases and three prostate cases were also used to demonstrate the advantages of our method. Our multi-GPU implementation can finish the optimization within ~1 minute for the H&N patient case. S1 leads to an inferior plan quality although its total time was 10 seconds shorter than the multi-GPU implementation. S2 and S3 yield same plan quality as the multi-GPU implementation but take ~4 minutes and ~6 minutes, respectively. High computational efficiency was consistently achieved for the other 5 cases. The results demonstrate that the multi-GPU implementation can handle the large-scale VMAT optimization problem efficiently without sacrificing plan quality., Comment: 23 pages, 3 figures, 3 tables

Monte Carlo (MC) simulation is considered as the most accurate method for radiation dose calculations. Accuracy of a source model for a linear accelerator is critical for the overall dose calculation accuracy. In this paper, we presented an analytical source model that we recently developed for GPU-based MC dose calculations. A key concept called phase-space-ring (PSR) was proposed. It contained a group of particles that are of the same type and close in energy and radial distance to the center of the phase-space plane. The model parameterized probability densities of particle location, direction and energy for each primary photon PSR, scattered photon PSR and electron PSR. For a primary photon PSRs, the particle direction is assumed to be from the beam spot. A finite spot size is modeled with a 2D Gaussian distribution. For a scattered photon PSR, multiple Gaussian components were used to model the particle direction. The direction distribution of an electron PSRs was also modeled as a 2D Gaussian distribution with a large standard deviation. We also developed a method to analyze a phase-space file and derive corresponding model parameters. To test the accuracy of our linac source model, dose distributions of different open fields in a water phantom were calculated using our source model and compared to those directly calculated using the reference phase-space file. The average distance-to-agreement (DTA) was within 1 mm for the depth dose in the build-up region and beam penumbra regions. The root-mean-square (RMS) dose difference was within 1.1% for dose profiles at inner and outer beam regions. The maximal relative difference of output factors was within 0.5%. Good agreements were also found in an IMRT prostate patient case and an IMRT head-and-neck case. These results demonstrated the efficacy of our source model in terms of accurately representing a reference phase-space file., Comment: page line has been removed

Monte Carlo (MC) method has been recognized the most accurate dose calculation method for radiotherapy. However, its extremely long computation time impedes clinical applications. Recently, a lot of efforts have been made to realize fast MC dose calculation on GPUs. Nonetheless, most of the GPU-based MC dose engines were developed in NVidia CUDA environment. This limits the code portability to other platforms, hindering the introduction of GPU-based MC simulations to clinical practice. The objective of this paper is to develop a fast cross-platform MC dose engine oclMC using OpenCL environment for external beam photon and electron radiotherapy in MeV energy range. Coupled photon-electron MC simulation was implemented with analogue simulations for photon transports and a Class II condensed history scheme for electron transports. To test the accuracy and efficiency of our dose engine oclMC, we compared dose calculation results of oclMC and gDPM, our previously developed GPU-based MC code, for a 15 MeV electron beam and a 6 MV photon beam on a homogenous water phantom, one slab phantom and one half-slab phantom. Satisfactory agreement was observed in all the cases. The average dose differences within 10% isodose line of the maximum dose were 0.48-0.53% for the electron beam cases and 0.15-0.17% for the photon beam cases. In terms of efficiency, our dose engine oclMC was 6-17% slower than gDPM when running both codes on the same NVidia TITAN card due to both different physics particle transport models and different computational environments between CUDA and OpenCL. The cross-platform portability was also validated by successfully running our new dose engine on a set of different compute devices including an Nvidia GPU card, two AMD GPU cards and an Intel CPU card using one or four cores. Computational efficiency among these platforms was compared., Comment: page line has been removed

In this paper, we present a new method to generate an instantaneous volumetric image using a single x-ray projection. To fully extract motion information hidden in projection images, we partitioned a projection image into small patches. We utilized a sparse learning method to automatically select patches that have a high correlation with principal component analysis (PCA) coefficients of a lung motion model. A model that maps the patch intensity to the PCA coefficients is built along with the patch selection process. Based on this model, a measured projection can be used to predict the PCA coefficients, which are further used to generate a motion vector field and hence a volumetric image. We have also proposed an intensity baseline correction method based on the partitioned projection, where the first and the second moments of pixel intensities at a patch in a simulated image are matched with those in a measured image via a linear transformation. The proposed method has been valid in simulated data and real phantom data. The algorithm is able to identify patches that contain relevant motion information, e.g. diaphragm region. It is found that intensity correction step is important to remove the systematic error in the motion prediction. For the simulation case, the sparse learning model reduced prediction error for the first PCA coefficient to 5%, compared to the 10% error when sparse learning is not used. 95th percentile error for the predicted motion vector is reduced from 2.40 mm to 0.92mm. In the phantom case, the predicted tumor motion trajectory is successfully reconstructed with 0.82 mm mean vector field error compared to 1.66 mm error without using the sparse learning method. The algorithm robustness with respect to sparse level, patch size, and existence of diaphragm, as well as computation time, has also been studied., Comment: 20 pages, 9 figures

Cone beam CT (CBCT) has been widely used for patient setup in image guided radiation therapy (IGRT). Radiation dose from CBCT scans has become a clinical concern. The purposes of this study are 1) to commission a GPU-based Monte Carlo (MC) dose calculation package gCTD for Varian On-Board Imaging (OBI) system and test the calculation accuracy, and 2) to quantitatively evaluate CBCT dose from the OBI system in typical IGRT scan protocols. We first conducted dose measurements in a water phantom. X-ray source model parameters used in gCTD are obtained through a commissioning process. gCTD accuracy is demonstrated by comparing calculations with measurements in water and in CTDI phantoms. 25 brain cancer patients are used to study dose in a standard-dose head protocol, and 25 prostate cancer patients are used to study dose in pelvis protocol and pelvis spotlight protocol. Mean dose to each organ is calculated. Mean dose to 2% voxels that have the highest dose is also computed to quantify the maximum dose. It is found that the mean dose value to an organ varies largely among patients. Moreover, dose distribution is highly non-homogeneous inside an organ. The maximum dose is found to be 1~3 times higher than the mean dose depending on the organ, and is up to 8 times higher for the entire body due to the very high dose region in bony structures. High computational efficiency has also been observed in our studies, such that MC dose calculation time is less than 5 min for a typical case., Comment: 19 pages, 6 figures, 7 tables

In the treatment plan optimization for intensity modulated radiation therapy (IMRT), dose-deposition coefficient (DDC) matrix is often pre-computed to parameterize the dose contribution to each voxel in the volume of interest from each beamlet of unit intensity. However, due to the limitation of computer memory and the requirement on computational efficiency, in practice matrix elements of small values are usually truncated, which inevitably compromises the quality of the resulting plan. A fixed-point iteration scheme has been applied in IMRT optimization to solve this problem, which has been reported to be effective and efficient based on the observations of the numerical experiments. In this paper, we aim to point out the mathematics behind this scheme and to answer the following three questions: 1) whether the fixed-point iteration algorithm converges or not? 2) when it converges, whether the fixed point solution is same as the original solution obtained with the complete DDC matrix? 3) if not the same, whether the fixed point solution is more accurate than the naive solution of the truncated problem obtained without the fixed-point iteration? To answer these questions, we first performed mathematical analysis and deductions using a simplified fluence map optimization (FMO) model. Then we conducted numerical experiments on a head-and-neck patient case using both the simplified and the original FMO model. Both our mathematical analysis and numerical experiments demonstrate that with proper DDC matrix truncation, the fixed-point iteration can converge. Even though the converged solution is not the one that we obtain with the complete DDC matrix, the fixed-point iteration scheme could significantly improve the plan accuracy compared with the solution to the truncated problem obtained without the fixed-point iteration., Comment: 16 pages, 3 figures, 2 tables

A novel phase-space source implementation has been designed for GPU-based Monte Carlo dose calculation engines. Due to the parallelized nature of GPU hardware, it is essential to simultaneously transport particles of the same type and similar energies but separated spatially to yield a high efficiency. We present three methods for phase-space implementation that have been integrated into the most recent version of the GPU-based Monte Carlo radiotherapy dose calculation package gDPM v3.0. The first method is to sequentially read particles from a patient-dependent phase-space and sort them on-the-fly based on particle type and energy. The second method supplements this with a simple secondary collimator model and fluence map implementation so that patient-independent phase-space sources can be used. Finally, as the third method (called the phase-space-let, or PSL, method) we introduce a novel strategy to pre-process patient-independent phase-spaces and bin particles by type, energy and position. Position bins located outside a rectangular region of interest enclosing the treatment field are ignored, substantially decreasing simulation time. The three methods were validated in absolute dose against BEAMnrc/DOSXYZnrc and compared using gamma-index tests (2%/2mm above the 10% isodose). It was found that the PSL method has the optimal balance between accuracy and efficiency and thus is used as the default method in gDPM v3.0. Using the PSL method, open fields of 4x4, 10x10 and 30x30 cm2 in water resulted in gamma passing rates of 99.96%, 99.92% and 98.66%, respectively. Relative output factors agreed within 1%. An IMRT patient plan using the PSL method resulted in a passing rate of 97%, and was calculated in 50 seconds using a single GPU compared to 8.4 hours (per CPU) for BEAMnrc/DOSXYZnrc., Comment: 19 pages, 6 figures

In adaptive radiotherapy, deformable image registration is often conducted between the planning CT and treatment CT (or cone beam CT) to generate a deformation vector field (DVF) for dose accumulation and contour propagation. The auto propagated contours on the treatment CT may contain relatively large errors, especially in low contrast regions. A clinician inspection and editing of the propagated contours are frequently needed. The edited contours are able to meet the clinical requirement for adaptive therapy; however, the DVF is still inaccurate and inconsistent with the edited contours. The purpose of this work is to develop a contour-guided deformable image registration (CG-DIR) algorithm to improve the accuracy and consistency of the DVF for adaptive radiotherapy. Incorporation of the edited contours into the registration algorithm is realized by regularizing the objective function of the original demons algorithm with a term of intensity matching between the delineated structures set pairs. The CG-DIR algorithm is implemented on computer graphics processing units (GPUs) by following the original GPU-based demons algorithm computation framework [Gu et al, Phys Med Biol. 55(1): 207-219, 2010]. The performance of CG-DIR is evaluated on five clinical head-and-neck and one pelvic cancer patient data. It is found that compared with the original demons, CG-DIR improves the accuracy and consistency of the DVF, while retaining similar high computational efficiency.

Patient respiratory signal associated with the cone beam CT (CBCT) projections is important for lung cancer radiotherapy. In contrast to monitoring an external surrogate of respiration, such signal can be extracted directly from the CBCT projections. In this paper, we propose a novel local principle component analysis (LPCA) method to extract the respiratory signal by distinguishing the respiration motion-induced content change from the gantry rotation-induced content change in the CBCT projections. The LPCA method is evaluated by comparing with three state-of-the-art projection-based methods, namely, the Amsterdam Shroud (AS) method, the intensity analysis (IA) method, and the Fourier-transform based phase analysis (FT-p) method. The clinical CBCT projection data of eight patients, acquired under various clinical scenarios, were used to investigate the performance of each method. We found that the proposed LPCA method has demonstrated the best overall performance for cases tested and thus is a promising technique for extracting respiratory signal. We also identified the applicability of each existing method., Comment: 21 pages, 11 figures, submitted to Phys. Med. Biol

In a treatment plan optimization problem for radiotherapy, a clinically acceptable plan is usually generated by an optimization process with weighting factors or reference doses adjusted for organs. Recent discoveries indicate that adjusting parameters associated with each voxel may lead to better plan quality. However, it is still unclear regarding the mathematical reasons behind it. To answer questions related to this problem, we establish in this work a new mathematical framework equipped with two theorems. The new framework clarifies the different consequences of adjusting organ-dependent and voxel-dependent parameters for the treatment plan optimization of radiation therapy, as well as the different effects of adjusting weighting factors versus reference doses in the optimization process. The main discoveries are threefold: 1) While in the organ-based model the selection of the objective function has an impact on the quality of the optimized plans, this is no longer an issue for the voxel-based model since the entire Pareto surface could be generated regardless the specific form of the objective function as long as it satisfies certain mathematical conditions; 2) A larger Pareto surface is explored by adjusting voxel-dependent parameters than by adjusting organ-dependent parameters, possibly allowing for the generation of plans with better trade-offs among different clinical objectives; 3) Adjusting voxel weighting factors is preferred to adjusting the voxel reference doses since the Pareto optimality can be maintained.

The gamma-index test has been commonly adopted to quantify the degree of agreement between a reference dose distribution and an evaluation dose distribution. Monte Carlo (MC) simulation has been widely used for the radiotherapy dose calculation for both clinical and research purposes. The goal of this work is to investigate both theoretically and experimentally the impact of the MC statistical fluctuation on the gamma-index test when the fluctuation exists in the reference, the evaluation, or both dose distributions. To the first order approximation, we theoretically demonstrated in a simplified model that the statistical fluctuation tends to overestimate gamma-index values when existing in the reference dose distribution and underestimate gamma-index values when existing in the evaluation dose distribution given the original gamma-index is relatively large for the statistical fluctuation. Our numerical experiments using clinical photon radiation therapy cases have shown that 1) when performing a gamma-index test between an MC reference dose and a non-MC evaluation dose, the average gamma-index is overestimated and the passing rate decreases with the increase of the noise level in the reference dose; 2) when performing a gamma-index test between a non-MC reference dose and an MC evaluation dose, the average gamma-index is underestimated when they are within the clinically relevant range and the passing rate increases with the increase of the noise level in the evaluation dose; 3) when performing a gamma-index test between an MC reference dose and an MC evaluation dose, the passing rate is overestimated due to the noise in the evaluation dose and underestimated due to the noise in the reference dose. We conclude that the gamma-index test should be used with caution when comparing dose distributions computed with Monte Carlo simulation.

Simulation of x-ray projection images plays an important role in cone beam CT (CBCT) related research projects. A projection image contains primary signal, scatter signal, and noise. It is computationally demanding to perform accurate and realistic computations for all of these components. In this work, we develop a package on GPU, called gDRR, for the accurate and efficient computations of x-ray projection images in CBCT under clinically realistic conditions. The primary signal is computed by a tri-linear ray-tracing algorithm. A Monte Carlo (MC) simulation is then performed, yielding the primary signal and the scatter signal, both with noise. A denoising process is applied to obtain a smooth scatter signal. The noise component is then obtained by combining the difference between the MC primary and the ray-tracing primary signals, and the difference between the MC simulated scatter and the denoised scatter signals. Finally, a calibration step converts the calculated noise signal into a realistic one by scaling its amplitude. For a typical CBCT projection with a poly-energetic spectrum, the calculation time for the primary signal is 1.2~2.3 sec, while the MC simulations take 28.1~95.3 sec. Computation time for all other steps is negligible. The ray-tracing primary signal matches well with the primary part of the MC simulation result. The MC simulated scatter signal using gDRR is in agreement with EGSnrc results with a relative difference of 3.8%. A noise calibration process is conducted to calibrate gDRR against a real CBCT scanner. The calculated projections are accurate and realistic, such that beam-hardening artifacts and scatter artifacts can be reproduced using the simulated projections. The noise amplitudes in the CBCT images reconstructed from the simulated projections also agree with those in the measured images at corresponding mAs levels., Comment: 21 pages, 11 figures, 1 table

Computed tomography (CT) to cone-beam computed tomography (CBCT) deformable image registration (DIR) is a crucial step in adaptive radiation therapy. Current intensity-based registration algorithms, such as demons, may fail in the context of CT-CBCT DIR because of inconsistent intensities between the two modalities. In this paper, we propose a variant of demons, called Deformation with Intensity Simultaneously Corrected (DISC), to deal with CT-CBCT DIR. DISC distinguishes itself from the original demons algorithm by performing an adaptive intensity correction step on the CBCT image at every iteration step of the demons registration. Specifically, the intensity correction of a voxel in CBCT is achieved by matching the first and the second moments of the voxel intensities inside a patch around the voxel with those on the CT image. It is expected that such a strategy can remove artifacts in the CBCT image, as well as ensuring the intensity consistency between the two modalities. DISC is implemented on computer graphics processing units (GPUs) in compute unified device architecture (CUDA) programming environment. The performance of DISC is evaluated on a simulated patient case and six clinical head-and-neck cancer patient data. It is found that DISC is robust against the CBCT artifacts and intensity inconsistency and significantly improves the registration accuracy when compared with the original demons., Comment: 22 pages, 9 figures, and 3 tables

Using fiducial markers on patient's body surface to predict the tumor location is a widely used approach in lung cancer radiotherapy. The purpose of this work is to propose an algorithm that automatically identifies a sparse set of locations on the patient's surface with the optimal prediction power for the tumor motion. The sparse selection of markers on the external surface and the assumed linear relationship between the marker motion and the internal tumor motion are represented by a prediction matrix. Such a matrix is determined by solving an optimization problem, where the objective function contains a sparsity term that penalizes the number of markers chosen on the patient's surface. The performance of our algorithm has been tested on realistic clinical data of four lung cancer patients. Thoracic 4DCT scans with 10 phases are used for the study. On a reference phase, a grid of points are casted on the patient's surface (except for patient's back) and propagated to other phases via deformable image registration of the corresponding CT images. Tumor locations at each phase are also manually delineated. We use 9 out of 10 phases of the 4DCT images to identify a small group of surface markers that are most correlated with the motion of the tumor, and find the prediction matrix at the same time. The 10th phase is then used to test the accuracy of the prediction. It is found that on average 6 to 7 surface markers are necessary to predict tumor locations with a 3D error of about 1mm. In addition, the selected marker locations lie closely in those areas where surface point motion has a high correlation with the tumor motion. Our method can automatically select sparse locations on patient's external surface and estimate a correlation matrix based on 4DCT, so that the selected surface locations can be used to place fiducial markers to optimally predict internal tumor motions.

Respiration-correlated CBCT, commonly called 4DCBCT, provide respiratory phase-resolved CBCT images. In many clinical applications, it is more preferable to reconstruct true 4DCBCT with the 4th dimension being time, i.e., each CBCT image is reconstructed based on the corresponding instantaneous projection. We propose in this work a novel algorithm for the reconstruction of this truly time-resolved CBCT, called cine-CBCT, by effectively utilizing the underlying temporal coherence, such as periodicity or repetition, in those cine-CBCT images. Assuming each column of the matrix $\bm{U}$ represents a CBCT image to be reconstructed and the total number of columns is the same as the number of projections, the central idea of our algorithm is that the rank of $\bm{U}$ is much smaller than the number of projections and we can use a matrix factorization form $\bm{U}=\bm{L}\bm{R}$ for $\bm{U}$. The number of columns for the matrix $\bm{L}$ constraints the rank of $\bm{U}$ and hence implicitly imposing a temporal coherence condition among all the images in cine-CBCT. The desired image properties in $\bm{L}$ and the periodicity of the breathing pattern are achieved by penalizing the sparsity of the tight wavelet frame transform of $\bm{L}$ and that of the Fourier transform of $\bm{R}$, respectively. A split Bregman method is used to solve the problem. In this paper we focus on presenting this new algorithm and showing the proof of principle using simulation studies on an NCAT phantom., Comment: 17 pages, 6 figures

Four-dimensional Cone Beam Computed Tomography (4D-CBCT) has been developed to provide respiratory phase resolved volumetric imaging in image guided radiation therapy (IGRT). Inadequate number of projections in each phase bin results in low quality 4D-CBCT images with obvious streaking artifacts. In this work, we propose two novel 4D-CBCT algorithms: an iterative reconstruction algorithm and an enhancement algorithm, utilizing a temporal nonlocal means (TNLM) method. We define a TNLM energy term for a given set of 4D-CBCT images. Minimization of this term favors those 4D-CBCT images such that any anatomical features at one spatial point at one phase can be found in a nearby spatial point at neighboring phases. 4D-CBCT reconstruction is achieved by minimizing a total energy containing a data fidelity term and the TNLM energy term. As for the image enhancement, 4D-CBCT images generated by the FDK algorithm are enhanced by minimizing the TNLM function while keeping the enhanced images close to the FDK results. A forward-backward splitting algorithm and a Gauss-Jacobi iteration method are employed to solve the problems. The algorithms are implemented on GPU to achieve a high computational efficiency. The reconstruction algorithm and the enhancement algorithm generate visually similar 4D-CBCT images, both better than the FDK results. Quantitative evaluations indicate that, compared with the FDK results, our reconstruction method improves contrast-to-noise-ratio (CNR) by a factor of 2.56~3.13 and our enhancement method increases the CNR by 2.75~3.33 times. The enhancement method also removes over 80% of the streak artifacts from the FDK results. The total computation time is ~460 sec for the reconstruction algorithm and ~610 sec for the enhancement algorithm on an NVIDIA Tesla C1060 GPU card., Comment: 20 pages, 3 figures, 2 tables

While compressed sensing (CS) based reconstructions have been developed for low-dose CBCT, a clear understanding on the relationship between the image quality and imaging dose at low dose levels is needed. In this paper, we qualitatively investigate this subject in a comprehensive manner with extensive experimental and simulation studies. The basic idea is to plot image quality and imaging dose together as functions of number of projections and mAs per projection over the whole clinically relevant range. A clear understanding on the tradeoff between image quality and dose can be achieved and optimal low-dose CBCT scan protocols can be developed for various imaging tasks in IGRT. Main findings of this work include: 1) Under the CS framework, image quality has little degradation over a large dose range, and the degradation becomes evident when the dose < 100 total mAs. A dose < 40 total mAs leads to a dramatic image degradation. Optimal low-dose CBCT scan protocols likely fall in the dose range of 40-100 total mAs, depending on the specific IGRT applications. 2) Among different scan protocols at a constant low-dose level, the super sparse-view reconstruction with projection number less than 50 is the most challenging case, even with strong regularization. Better image quality can be acquired with other low mAs protocols. 3) The optimal scan protocol is the combination of a medium number of projections and a medium level of mAs/view. This is more evident when the dose is around 72.8 total mAs or below and when the ROI is a low-contrast or high-resolution object. Based on our results, the optimal number of projections is around 90 to 120. 4) The clinically acceptable lowest dose level is task dependent. In our study, 72.8mAs is a safe dose level for visualizing low-contrast objects, while 12.2 total mAs is sufficient for detecting high-contrast objects of diameter greater than 3 mm., Comment: 19 pages, 12 figures, submitted to Physics in Medicine and Biology

Recently, X-ray imaging dose from computed tomography (CT) or cone beam CT (CBCT) scans has become a serious concern. Patient-specific imaging dose calculation has been proposed for the purpose of dose management. While Monte Carlo (MC) dose calculation can be quite accurate for this purpose, it suffers from low computational efficiency. In response to this problem, we have successfully developed a MC dose calculation package, gCTD, on GPU architecture under the NVIDIA CUDA platform for fast and accurate estimation of the x-ray imaging dose received by a patient during a CT or CBCT scan. Techniques have been developed particularly for the GPU architecture to achieve high computational efficiency. Dose calculations using CBCT scanning geometry in a homogeneous water phantom and a heterogeneous Zubal head phantom have shown good agreement between gCTD and EGSnrc, indicating the accuracy of our code. In terms of improved efficiency, it is found that gCTD attains a speed-up of ~400 times in the homogeneous water phantom and ~76.6 times in the Zubal phantom compared to EGSnrc. As for absolute computation time, imaging dose calculation for the Zubal phantom can be accomplished in ~17 sec with the average relative standard deviation of 0.4%. Though our gCTD code has been developed and tested in the context of CBCT scans, with simple modification of geometry it can be used for assessing imaging dose in CT scans as well., Comment: 18 pages, 7 figures, and 1 table

Monte Carlo (MC) simulation is commonly considered to be the most accurate dose calculation method in radiotherapy. However, its efficiency still requires improvement for many routine clinical applications. In this paper, we present our recent progress towards the development a GPU-based MC dose calculation package, gDPM v2.0. It utilizes the parallel computation ability of a GPU to achieve high efficiency, while maintaining the same particle transport physics as in the original DPM code and hence the same level of simulation accuracy. In GPU computing, divergence of execution paths between threads can considerably reduce the efficiency. Since photons and electrons undergo different physics and hence attain different execution paths, we use a simulation scheme where photon transport and electron transport are separated to partially relieve the thread divergence issue. High performance random number generator and hardware linear interpolation are also utilized. We have also developed various components to handle fluence map and linac geometry, so that gDPM can be used to compute dose distributions for realistic IMRT or VMAT treatment plans. Our gDPM package is tested for its accuracy and efficiency in both phantoms and realistic patient cases. In all cases, the average relative uncertainties are less than 1%. A statistical t-test is performed and the dose difference between the CPU and the GPU results is found not statistically significant in over 96% of the high dose region and over 97% of the entire region. Speed up factors of 69.1 ~ 87.2 have been observed using an NVIDIA Tesla C2050 GPU card against a 2.27GHz Intel Xeon CPU processor. For realistic IMRT and VMAT plans, MC dose calculation can be completed with less than 1% standard deviation in 36.1~39.6 sec using gDPM., Comment: 18 pages, 5 figures, and 3 tables

Targeting at the development of an accurate and efficient dose calculation engine for online adaptive radiotherapy, we have implemented a finite size pencil beam (FSPB) algorithm with a 3D-density correction method on GPU. This new GPU-based dose engine is built on our previously published ultrafast FSPB computational framework [Gu et al. Phys. Med. Biol. 54 6287-97, 2009]. Dosimetric evaluations against Monte Carlo dose calculations are conducted on 10 IMRT treatment plans (5 head-and-neck cases and 5 lung cases). For all cases, there is improvement with the 3D-density correction over the conventional FSPB algorithm and for most cases the improvement is significant. Regarding the efficiency, because of the appropriate arrangement of memory access and the usage of GPU intrinsic functions, the dose calculation for an IMRT plan can be accomplished well within 1 second (except for one case) with this new GPU-based FSPB algorithm. Compared to the previous GPU-based FSPB algorithm without 3D-density correction, this new algorithm, though slightly sacrificing the computational efficiency (~5-15% lower), has significantly improved the dose calculation accuracy, making it more suitable for online IMRT replanning.

Recently we have developed an algorithm for reconstructing volumetric images and extracting 3D tumor motion information from a single x-ray projection. We have demonstrated its feasibility using a digital respiratory phantom with regular breathing patterns. In this work, we present a detailed description and a comprehensive evaluation of the improved algorithm. The algorithm was improved by incorporating respiratory motion prediction. The accuracy and efficiency were then evaluated on 1) a digital respiratory phantom, 2) a physical respiratory phantom, and 3) five lung cancer patients. These evaluation cases include both regular and irregular breathing patterns that are different from the training dataset. For the digital respiratory phantom with regular and irregular breathing, the average 3D tumor localization error is less than 1 mm. On an NVIDIA Tesla C1060 GPU card, the average computation time for 3D tumor localization from each projection ranges between 0.19 and 0.26 seconds, for both regular and irregular breathing, which is about a 10% improvement over previously reported results. For the physical respiratory phantom, an average tumor localization error below 1 mm was achieved with an average computation time of 0.13 and 0.16 seconds on the same GPU card, for regular and irregular breathing, respectively. For the five lung cancer patients, the average tumor localization error is below 2 mm in both the axial and tangential directions. The average computation time on the same GPU card ranges between 0.26 and 0.34 seconds.

Beam orientation optimization (BOO) is a key component in the process of IMRT treatment planning. It determines to what degree one can achieve a good treatment plan quality in the subsequent plan optimization process. In this paper, we have developed a BOO algorithm via adaptive l_1 minimization. Specifically, we introduce a sparsity energy function term into our model which contains weighting factors for each beam angle adaptively adjusted during the optimization process. Such an energy term favors small number of beam angles. By optimizing a total energy function containing a dosimetric term and the sparsity term, we are able to identify the unimportant beam angles and gradually remove them without largely sacrificing the dosimetric objective. In one typical prostate case, the convergence property of our algorithm, as well as the how the beam angles are selected during the optimization process, is demonstrated. Fluence map optimization (FMO) is then performed based on the optimized beam angles. The resulted plan quality is presented and found to be better than that obtained from unoptimized (equiangular) beam orientations. We have further systematically validated our algorithm in the contexts of 5-9 coplanar beams for 5 prostate cases and 1 head and neck case. For each case, the final FMO objective function value is used to compare the optimized beam orientations and the equiangular ones. It is found that, our BOO algorithm can lead to beam configurations which attain lower FMO objective function values than corresponding equiangular cases, indicating the effectiveness of our BOO algorithm., Comment: 19 pages, 2 tables, and 5 figures

The gamma-index dose comparison tool has been widely used to compare dose distributions in cancer radiotherapy. The accurate calculation of gamma-index requires an exhaustive search of the closest Euclidean distance in the high-resolution dose-distance space. This is a computational intensive task when dealing with 3D dose distributions. In this work, we combine a geometric method with a radial pre-sorting technique , and implement them on computer graphics processing units (GPUs). The developed GPU-based gamma-index computational tool is evaluated on eight pairs of IMRT dose distributions. The GPU implementation achieved 20x~30x speedup factor compared to CPU implementation and gamma-index calculations can be finished within a few seconds for all 3D testing cases. We further investigated the effect of various factors on both CPU and GPU computation time. The strategy of pre-sorting voxels based on their dose difference values speed up the GPU calculation by about 2-4 times. For n-dimensional dose distributions, gamma-index calculation time on CPU is proportional to the summation of gamma^n over all voxels, while that on GPU is effected by gamma^n distributions and is approximately proportional to the gamma^n summation over all voxels. We found increasing dose distributions resolution leads to quadratic increase of computation time on CPU, while less-than-quadratic increase on GPU. The values of dose difference (DD) and distance-to-agreement (DTA) criteria also have their impact on gamma-index calculation time., Comment: 13 pages, 2 figures, and 3 tables

High radiation dose in CT scans increases a lifetime risk of cancer and has become a major clinical concern. Recently, iterative reconstruction algorithms with Total Variation (TV) regularization have been developed to reconstruct CT images from highly undersampled data acquired at low mAs levels in order to reduce the imaging dose. Nonetheless, TV regularization may lead to over-smoothed images and lost edge information. To solve this problem, in this work we develop an iterative CT reconstruction algorithm with edge-preserving TV regularization to reconstruct CT images from highly undersampled data obtained at low mAs levels. The CT image is reconstructed by minimizing an energy consisting of an edge-preserving TV norm and a data fidelity term posed by the x-ray projections. The edge-preserving TV term is proposed to preferentially perform smoothing only on non-edge part of the image in order to avoid over-smoothing, which is realized by introducing a penalty weight to the original total variation norm. Our iterative algorithm is implemented on GPU to improve its speed. We test our reconstruction algorithm on a digital NCAT phantom, a physical chest phantom, and a Catphan phantom. Reconstruction results from a conventional FBP algorithm and a TV regularization method without edge preserving penalty are also presented for comparison purpose. The experimental results illustrate that both TV-based algorithm and our edge-preserving TV algorithm outperform the conventional FBP algorithm in suppressing the streaking artifacts and image noise under the low dose context. Our edge-preserving algorithm is superior to the TV-based algorithm in that it can preserve more information of fine structures and therefore maintain acceptable spatial resolution., Comment: 21 pages, 6 figures, 2 tables

Four-dimensional computed tomography (4DCT) has been widely used in cancer radiotherapy for accurate target delineation and motion measurement for tumors in thorax and upper abdomen areas. However, 4DCT simulation is associated with much higher imaging dose than conventional CT simulation, which is a major concern in its clinical application. Conventionally, each phase of 4DCT is reconstructed independently using the filtered backprojection (FBP) algorithm. The basic idea of our new algorithm is that, by utilizing the common information among different phases, the input information required to reconstruct image of high quality, and thus the imaging dose, can be reduced. We proposed a temporal non-local means (TNLM) method to explore the inter-phase similarity. All phases of the 4DCT images are reconstructed simultaneously by minimizing a cost function consisting of a data fidelity term and a TNLM regularization term. We utilized a forward-backward splitting algorithm and a Gauss-Jacobi iteration method to efficiently solve the minimization problem. The algorithm was also implemented on graphics processing unit (GPU) to achieve a high computational speed. Our reconstruction algorithm has been tested on a digital NCAT thorax phantom in three low dose scenarios. Our new algorithm generates visually much better CT images containing less image noise and streaking artifacts compared with the standard FBP algorithm. Quantitative analysis shows that much higher contrast-to-noise ratio and signal-to-noise ratio can be obtained using our algorithm. The total reconstruction time for all 10 phases of a slice ranges from 90 to 140 seconds on an NVIDIA Tesla C1060 GPU card., Comment: 14 pages, 4 figures, and 2 tables, accepted by Med. Phys

X-ray imaging dose from serial cone-beam CT (CBCT) scans raises a clinical concern in most image guided radiation therapy procedures. It is the goal of this paper to develop a fast GPU-based algorithm to reconstruct high quality CBCT images from undersampled and noisy projection data so as to lower the imaging dose. For this purpose, we have developed an iterative tight frame (TF) based CBCT reconstruction algorithm. A condition that a real CBCT image has a sparse representation under a TF basis is imposed in the iteration process as regularization to the solution. To speed up the computation, a multi-grid method is employed. Our GPU implementation has achieved high computational efficiency and a CBCT image of resolution 512\times512\times70 can be reconstructed in ~5 min. We have tested our algorithm on a digital NCAT phantom and a physical Catphan phantom. It is found that our TF-based algorithm is able to reconstrct CBCT in the context of undersampling and low mAs levels. We have also quantitatively analyzed the reconstructed CBCT image quality in terms of modulation-transfer-function and contrast-to-noise ratio under various scanning conditions. The results confirm the high CBCT image quality obtained from our TF algorithm. Moreover, our algorithm has also been validated in a real clinical context using a head-and-neck patient case. Comparisons of the developed TF algorithm and the current state-of-the-art TV algorithm have also been made in various cases studied in terms of reconstructed image quality and computation efficiency., Comment: 24 pages, 8 figures, accepted by Phys. Med. Biol

Purpose: To develop a novel aperture-based algorithm for volumetric modulated arc therapy (VMAT) treatment plan optimization with high quality and high efficiency. Methods: The VMAT optimization problem is formulated as a large-scale convex programming problem solved by a column generation approach. We consider a cost function consisting two terms, the first which enforces a desired dose distribution while the second guarantees a smooth dose rate variation between successive gantry angles. At each iteration of the column generation method, a subproblem is first solved to generate one more deliverable MLC aperture which potentially decreases the cost function most effectively. A subsequent master problem is then solved to determine the dose rate at all currently generated apertures by minimizing the cost function. The iteration of such an algorithm yields a set of deliverable apertures, as well as dose rates, at all gantry angles. Results: The algorithm was preliminarily tested on five prostate and five head-and-neck clinical cases, each with one full gantry rotation and without any couch/collimator rotations. Compared to corresponding co-planar IMRT treatment plans (9 fields for prostate cases and 5 fields for head-and-neck cases), the VMAT plans delivered much lower doses to critical structures and more conformal doses to targets. Moreover, extremely high efficiency has been achieved in our algorithm. It takes only 5~8 minutes on CPU (MATLAB code on an Intel Xeon 2.27 GHz CPU) and 18~31 seconds on GPU (CUDA code on an NVIDIA Tesla C1060 GPU card) to generate such a plan. Conclusions: We have developed an aperture-based VMAT optimization algorithm which can generate clinically deliverable and high quality treatment plans at very high efficiency.

Purpose: To develop an algorithm for real-time volumetric image reconstruction and 3D tumor localization based on a single x-ray projection image for lung cancer radiotherapy. Methods: Given a set of volumetric images of a patient at N breathing phases as the training data, we perform deformable image registration between a reference phase and the other N-1 phases, resulting in N-1 deformation vector fields (DVFs). These DVFs can be represented efficiently by a few eigenvectors and coefficients obtained from principal component analysis (PCA). By varying the PCA coefficients, we can generate new DVFs, which, when applied on the reference image, lead to new volumetric images. We then can reconstruct a volumetric image from a single projection image by optimizing the PCA coefficients such that its computed projection matches the measured one. The 3D location of the tumor can be derived by applying the inverted DVF on its position in the reference image. Our algorithm was implemented on graphics processing units (GPUs) to achieve real-time efficiency. We generated the training data using a realistic and dynamic mathematical phantom with 10 breathing phases. The testing data were 360 cone beam projections corresponding to one gantry rotation, simulated using the same phantom with a 50% increase in breathing amplitude. Results: The average relative image intensity error of the reconstructed volumetric images is 6.9% +/- 2.4%. The average 3D tumor localization error is 0.8 mm +/- 0.5 mm. On an NVIDIA Tesla C1060 GPU card, the average computation time for reconstructing a volumetric image from each projection is 0.24 seconds (range: 0.17 and 0.35 seconds). Conclusions: We have shown the feasibility of reconstructing volumetric images and localizing tumor positions in 3D in near real-time from a single x-ray image., Comment: 8 pages, 3 figures, submitted to Medical Physics Letter

Online adaptive radiation therapy (ART) has great promise to significantly reduce normal tissue toxicity and/or improve tumor control through real-time treatment adaptations based on the current patient anatomy. However, the major technical obstacle for clinical realization of online ART, namely the inability to achieve real-time efficiency in treatment re-planning, has yet to be solved. To overcome this challenge, this paper presents our work on the implementation of an intensity modulated radiation therapy (IMRT) direct aperture optimization (DAO) algorithm on graphics processing unit (GPU) based on our previous work on CPU. We formulate the DAO problem as a large-scale convex programming problem, and use an exact method called column generation approach to deal with its extremely large dimensionality on GPU. Five 9-field prostate and five 5-field head-and-neck IMRT clinical cases with 5\times5 mm2 beamlet size and 2.5\times2.5\times2.5 mm3 voxel size were used to evaluate our algorithm on GPU. It takes only 0.7~2.5 seconds for our implementation to generate optimal treatment plans using 50 MLC apertures on an NVIDIA Tesla C1060 GPU card. Our work has therefore solved a major problem in developing ultra-fast (re-)planning technologies for online ART.

Purpose: Cone-beam CT (CBCT) plays an important role in image guided radiation therapy (IGRT). However, the large radiation dose from serial CBCT scans in most IGRT procedures raises a clinical concern, especially for pediatric patients who are essentially excluded from receiving IGRT for this reason. The goal of this work is to develop a fast GPU-based algorithm to reconstruct CBCT from undersampled and noisy projection data so as to lower the imaging dose. Methods: The CBCT is reconstructed by minimizing an energy functional consisting of a data fidelity term and a total variation regularization term. We developed a GPU-friendly version of the forward-backward splitting algorithm to solve this model. A multi-grid technique is also employed. Results: It is found that 20~40 x-ray projections are sufficient to reconstruct images with satisfactory quality for IGRT. The reconstruction time ranges from 77 to 130 sec on a NVIDIA Tesla C1060 GPU card, depending on the number of projections used, which is estimated about 100 times faster than similar iterative reconstruction approaches. Moreover, phantom studies indicate that our algorithm enables the CBCT to be reconstructed under a scanning protocol with as low as 0.1 mAs/projection. Comparing with currently widely used full-fan head and neck scanning protocol of ~360 projections with 0.4 mAs/projection, it is estimated that an overall 36~72 times dose reduction has been achieved in our fast CBCT reconstruction algorithm. Conclusions: This work indicates that the developed GPU-based CBCT reconstruction algorithm is capable of lowering imaging dose considerably. The high computation efficiency in this algorithm makes the iterative CBCT reconstruction approach applicable in real clinical environments., Comment: Accepted as a letter in Med. Phys., brief clarifying comments and updated references. 6 pages and 2 figures

Organ motion induced by respiration may cause clinically significant targeting errors and greatly degrade the effectiveness of conformal radiotherapy. It is therefore crucial to be able to model respiratory motion accurately. A recently proposed lung motion model based on principal component analysis (PCA) has been shown to be promising on a few patients. However, there is still a need to understand the underlying reason why it works. In this paper, we present a much deeper and detailed analysis of the PCA-based lung motion model. We provide the theoretical justification of the effectiveness of PCA in modeling lung motion. We also prove that under certain conditions, the PCA motion model is equivalent to 5D motion model, which is based on physiology and anatomy of the lung. The modeling power of PCA model was tested on clinical data and the average 3D error was found to be below 1 mm., Comment: 4 pages, 1 figure. submitted to International Conference on the use of Computers in Radiation Therapy 2010

Cone-beam CT (CBCT) has been widely used in image guided radiation therapy (IGRT) to acquire updated volumetric anatomical information before treatment fractions for accurate patient alignment purpose. However, the excessive x-ray imaging dose from serial CBCT scans raises a clinical concern in most IGRT procedures. The excessive imaging dose can be effectively reduced by reducing the number of x-ray projections and/or lowering mAs levels in a CBCT scan. The goal of this work is to develop a fast GPU-based algorithm to reconstruct high quality CBCT images from undersampled and noisy projection data so as to lower the imaging dose. The CBCT is reconstructed by minimizing an energy functional consisting of a data fidelity term and a total variation regularization term. We developed a GPU-friendly version of the forward-backward splitting algorithm to solve this model. A multi-grid technique is also employed. We test our CBCT reconstruction algorithm on a digital NCAT phantom and a head-and-neck patient case. The performance under low mAs is also validated using a physical Catphan phantom and a head-and-neck Rando phantom. It is found that 40 x-ray projections are sufficient to reconstruct CBCT images with satisfactory quality for IGRT patient alignment purpose. Phantom experiments indicated that CBCT images can be successfully reconstructed with our algorithm under as low as 0.1 mAs/projection level. Comparing with currently widely used full-fan head-and-neck scanning protocol of about 360 projections with 0.4 mAs/projection, it is estimated that an overall 36 times dose reduction has been achieved with our algorithm. Moreover, the reconstruction time is about 130 sec on an NVIDIA Tesla C1060 GPU card, which is estimated ~100 times faster than similar iterative reconstruction approaches., Comment: 20 pages, 10 figures, Paper was revised and more testing cases were added.

Monte Carlo simulation is the most accurate method for absorbed dose calculations in radiotherapy. Its efficiency still requires improvement for routine clinical applications, especially for online adaptive radiotherapy. In this paper, we report our recent development on a GPU-based Monte Carlo dose calculation code for coupled electron-photon transport. We have implemented the Dose Planning Method (DPM) Monte Carlo dose calculation package (Sempau et al, Phys. Med. Biol., 45(2000)2263-2291) on GPU architecture under CUDA platform. The implementation has been tested with respect to the original sequential DPM code on CPU in phantoms with water-lung-water or water-bone-water slab geometry. A 20 MeV mono-energetic electron point source or a 6 MV photon point source is used in our validation. The results demonstrate adequate accuracy of our GPU implementation for both electron and photon beams in radiotherapy energy range. Speed up factors of about 5.0 ~ 6.6 times have been observed, using an NVIDIA Tesla C1060 GPU card against a 2.27GHz Intel Xeon CPU processor., Comment: 13 pages, 3 figures, and 1 table. Paper revised. Figures updated

Online adaptive radiation therapy (ART) promises the ability to deliver an optimal treatment in response to daily patient anatomic variation. A major technical barrier for the clinical implementation of online ART is the requirement of rapid image segmentation. Deformable image registration (DIR) has been used as an automated segmentation method to transfer tumor/organ contours from the planning image to daily images. However, the current computational time of DIR is insufficient for online ART. In this work, this issue is addressed by using computer graphics processing units (GPUs). A grey-scale based DIR algorithm called demons and five of its variants were implemented on GPUs using the Compute Unified Device Architecture (CUDA) programming environment. The spatial accuracy of these algorithms was evaluated over five sets of pulmonary 4DCT images with an average size of 256x256x100 and more than 1,100 expert-determined landmark point pairs each. For all the testing scenarios presented in this paper, the GPU-based DIR computation required around 7 to 11 seconds to yield an average 3D error ranging from 1.5 to 1.8 mm. It is interesting to find out that the original passive force demons algorithms outperform subsequently proposed variants based on the combination of accuracy, efficiency, and ease of implementation., Comment: Submitted to Physics in Medicine and Biology

The widespread adoption of on-board volumetric imaging in cancer radiotherapy has stimulated research efforts to develop online adaptive radiotherapy techniques to handle the inter-fraction variation of the patient's geometry. Such efforts face major technical challenges to perform treatment planning in real-time. To overcome this challenge, we are developing a supercomputing online re-planning environment (SCORE) at the University of California San Diego (UCSD). As part of the SCORE project, this paper presents our work on the implementation of an intensity modulated radiation therapy (IMRT) optimization algorithm on graphics processing units (GPUs). We adopt a penalty-based quadratic optimization model, which is solved by using a gradient projection method with Armijo's line search rule. Our optimization algorithm has been implemented in CUDA for parallel GPU computing as well as in C for serial CPU computing for comparison purpose. A prostate IMRT case with various beamlet and voxel sizes was used to evaluate our implementation. On an NVIDIA Tesla C1060 GPU card, we have achieved speedup factors of 20-40 without losing accuracy, compared to the results from an Intel Xeon 2.27 GHz CPU. For a specific 9-field prostate IMRT case with 5x5 mm^2 beamlet size and 2.5x2.5x2.5 mm^3 voxel size, our GPU implementation takes only 2.8 seconds to generate an optimal IMRT plan. Our work has therefore solved a major problem in developing online re-planning technologies for adaptive radiotherapy., Comment: submitted to Physics in Medicine and Biology

Online adaptive radiation therapy (ART) is an attractive concept that promises the ability to deliver an optimal treatment in response to the inter-fraction variability in patient anatomy. However, it has yet to be realized due to technical limitations. Fast dose deposit coefficient calculation is a critical component of the online planning process that is required for plan optimization of intensity modulated radiation therapy (IMRT). Computer graphics processing units (GPUs) are well-suited to provide the requisite fast performance for the data-parallel nature of dose calculation. In this work, we develop a dose calculation engine based on a finite-size pencil beam (FSPB) algorithm and a GPU parallel computing framework. The developed framework can accommodate any FSPB model. We test our implementation on a case of a water phantom and a case of a prostate cancer patient with varying beamlet and voxel sizes. All testing scenarios achieved speedup ranging from 200~400 times when using a NVIDIA Tesla C1060 card in comparison with a 2.27GHz Intel Xeon CPU. The computational time for calculating dose deposition coefficients for a 9-field prostate IMRT plan with this new framework is less than 1 second. This indicates that the GPU-based FSPB algorithm is well-suited for online re-planning for adaptive radiotherapy., Comment: submitted Physics in Medicine and Biology

Lung tumor motion due to respiration poses a challenge in the application of modern three-dimensional conformal radiotherapy. Direct tracking of the lung tumor during radiation therapy is very difficult without implanted fiducial markers. Indirect tracking relies on the correlation of the tumor's motion and the surrogate's motion. The present paper presents an analysis of the correlation between the tumor motion and the diaphragm motion in order to evaluate the potential use of diaphragm as a surrogate for tumor motion. We have analyzed the correlation between diaphragm motion and superior-inferior lung tumor motion in 32 fluoroscopic image sequences from 10 lung cancer patients. A simple linear model and a more complex linear model that accounts for phase delays between the two motions have been used. Results show that the diaphragm is a good surrogate for tumor motion prediction for most patients, resulting in an average correlation factor of 0.94 and 0.98 with each model respectively. The model that accounts for delays leads to an average localization prediction error of 0.8mm and an error at the 95% confidence level of 2.1mm. However, for one patient studied, the correlation is much weaker compared to other patients. This indicates that, before using diaphragm for lung tumor prediction, the correlation should be examined on a patient-by-patient basis., Comment: Accepted by Physics in Medicine and Biology