Mingli Chen

A simple fixed-point approach to invert a deformation fielda)

Inversion of deformation fields is applied frequently to map images, dose, and contours between the reference frame and the study frame. A prevailing approach that takes the negative of the forward deformation as the inverse deformation is oversimplified and can cause large errors for large deformations or deformations that are composites of several deformations. Other approaches, including Newtons method and scatter data interpolation, either require the first derivative or are very inefficient. Here we propose an iterative approach that is easy to implement, converges quickly to the inverse when it does, and works for a majority of cases in practice. Our approach is rooted in fixed-point theory. We build a sequence to approximate the inverse deformation through iterative evaluation of the forward deformation. A sufficient but not necessary convergence condition (Lipschitz condition) and its proof are also given. Though this condition guarantees the convergence, it may not be met for an arbitrary deformation field. One should always check whether the inverse exists for the given forward deformation field by calculating its Jacobian. If nonpositive values of the Jacobian occur only for few voxels, this method will usually converge to a pseudoinverse. In case the iteration fails to converge, one should switch to other means of finding the inverse. We tested the proposed method on simulated 2D data and real 3D computed tomography data of a lung patient and compared our method with two implementations in the Insight Segmentation and Registration Toolkit (ITK). Typically less than ten iterations are needed for our method to get an inverse deformation field with clinically relevant accuracy. Based on the test results, our method is about ten times faster and yet ten times more accurate than ITKs iterative method for the same number of iterations. Simulations and real data tests demonstrated the efficacy and the accuracy of the proposed algorithm.

Ultrafast convolution/superposition using tabulated and exponential kernels on GPU

PURPOSE Collapsed-cone convolution/superposition (CCCS) dose calculation is the workhorse for IMRT dose calculation. The authors present a novel algorithm for computing CCCS dose on the modern graphic processing unit (GPU). METHODS The GPU algorithm includes a novel TERMA calculation that has no write-conflicts and has linear computation complexity. The CCCS algorithm uses either tabulated or exponential cumulative-cumulative kernels (CCKs) as reported in literature. The authors have demonstrated that the use of exponential kernels can reduce the computation complexity by order of a dimension and achieve excellent accuracy. Special attentions are paid to the unique architecture of GPU, especially the memory accessing pattern, which increases performance by more than tenfold. RESULTS As a result, the tabulated kernel implementation in GPU is two to three times faster than other GPU implementations reported in literature. The implementation of CCCS showed significant speedup on GPU over single core CPU. On tabulated CCK, speedups as high as 70 are observed; on exponential CCK, speedups as high as 90 are observed. CONCLUSIONS Overall, the GPU algorithm using exponential CCK is 1000-3000 times faster over a highly optimized single-threaded CPU implementation using tabulated CCK, while the dose differences are within 0.5% and 0.5 mm. This ultrafast CCCS algorithm will allow many time-sensitive applications to use accurate dose calculation.

Efficient gamma index calculation using fast Euclidean distance transform

The gamma index is a tool for dose distribution comparison. It combines both dose difference (DD) and distance to agreement (DTA) into a single quantity. Though it is an effective measure, making up for the inadequacy of DD or DTA alone, its calculation can be very time-consuming. For a k-D space with N quantization levels in each dimension, the complexity of the exhaustive search is O(N(2k)). In this work, we proposed an efficient method that reduces the complexity from O(N(2k)) to O(N(k)M), where M is the number of discretized dose values and is comparable to N. More precisely, by embedding the reference dose distribution in a (k+1)-D spatial-dose space, we can use fast Euclidean distance transform with linear complexity to obtain a table of gamma indices evaluated over a range of the (k+1)-D spatial-dose space. Then, to obtain gamma indices for the test dose distribution, it requires only table lookup with complexity O(N(k)). Such a table can also be used for other test dose distributions as long as the reference dose distribution is the same. Simulations demonstrated the efficiency of our proposed method. The speedup for 3D gamma index calculation is expected to be on the order of tens of thousands (from O(N(6)) to O(N(3)M)) if N is a few hundreds, which makes clinical usage of the 3D gamma index feasible. A byproduct of the gamma index table is that the gradient of the gamma index with respect to either the spatial or dose dimension can be easily derived. The gradient can be used to identify the main causes of the discrepancy from the reference distribution at any dose point in the test distribution or incorporated in treatment planning and machine parameter optimization.

Validation of GPU based TomoTherapy dose calculation engine.

PURPOSE The graphic processing unit (GPU) based TomoTherapy convolution/superposition(C/S) dose engine (GPU dose engine) achieves a dramatic performance improvement over the traditional CPU-cluster based TomoTherapy dose engine (CPU dose engine). Besides the architecture difference between the GPU and CPU, there are several algorithm changes from the CPU dose engine to the GPU dose engine. These changes made the GPU dose slightly different from the CPU-cluster dose. In order for the commercial release of the GPU dose engine, its accuracy has to be validated. METHODS Thirty eight TomoTherapy phantom plans and 19 patient plans were calculated with both dose engines to evaluate the equivalency between the two dose engines. Gamma indices (Γ) were used for the equivalency evaluation. The GPU dose was further verified with the absolute point dose measurement with ion chamber and film measurements for phantom plans. Monte Carlo calculation was used as a reference for both dose engines in the accuracy evaluation in heterogeneous phantom and actual patients. RESULTS The GPU dose engine showed excellent agreement with the current CPU dose engine. The majority of cases had over 99.99% of voxels with Γ(1%, 1 mm) < 1. The worst case observed in the phantom had 0.22% voxels violating the criterion. In patient cases, the worst percentage of voxels violating the criterion was 0.57%. For absolute point dose verification, all cases agreed with measurement to within ±3% with average error magnitude within 1%. All cases passed the acceptance criterion that more than 95% of the pixels have Γ(3%, 3 mm) < 1 in film measurement, and the average passing pixel percentage is 98.5%-99%. The GPU dose engine also showed similar degree of accuracy in heterogeneous media as the current TomoTherapy dose engine. CONCLUSIONS It is verified and validated that the ultrafast TomoTherapy GPU dose engine can safely replace the existing TomoTherapy cluster based dose engine without degradation in dose accuracy.

Dynamic tomotherapy delivery

PURPOSE Several dynamic techniques are introduced to speed up TomoTherapy delivery and improve longitudinal conformity. These techniques include dynamic jaws, dynamic couch, and their combinations. METHODS In general, dynamic jaws techniques allow jaws to move during a treatment. On the one hand, the jaws open wide to increase efficiency and thus reduce beam-on time. On the other hand, the jaws can close and follow the target border when sharp penumbra is required near the superior and inferior borders of tumor sites, which results in improved longitudinal dose conformity. The main purpose of the dynamic couch technique is to move the couch as fast as possible at variable speed to reduce beam-on time. Delivering most conformal dose as fast as possible requires a combination of dynamic jaws and dynamic couch techniques (DJDC). Motions of the jaws and couch are determined from the longitudinal fluence profile, which is calculated from an optimized leaf sinogram of small jaw width regular delivery or running start and stop delivery (RSS). We focused on RSS and DJDC in this study and also discussed other delivery techniques. RESULTS Several conceptual cases are simulated to compare different delivery techniques. The results show that beam-on time can be reduced by about 60% compared to regular delivery with a 2.5 cm jaw width (REG 2.5 cm) for these cases and arbitrary longitudinal fluence profiles can be delivered. Two clinical cases, a prostate and a head-and-neck case, with different delivery techniques are calculated. The results show that plan quality yielded by DJDC with a maximum 5.0 cm jaw width is overall comparable to or better than that of the existing REG 2.5 cm. CONCLUSIONS The DVH comparisons show better critical structure avoidance with the dynamic techniques. At the same time, beam-on time is reduced by about one half compared to REG 2.5 cm. Dynamic delivery techniques provide users more tools to speed up delivery and/or improve plan quality.

Theoretical analysis of the thread effect in helical TomoTherapy

PURPOSE The longitudinal dose ripple on the off-axis caused by helical radiation delivery, such as the TomoTherapy system, has been observed, and its relation with respect to pitch has been studied with empirically found optimal pitches, 0.86∕n, by Kissick et al. [Med. Phys. 32, 1414-1423 (2005)]. This ripple artifact referred to as the thread effect is periodic in nature and is caused by various periodic factors. In this work, the factors that cause the thread effect were unveiled, including jaw profile divergence, the inverse square law, attenuation, and the cone effect, and their impact on the thread effect were studied. METHODS Mathematical formulation for individual and combined factors were set up. Based on theoretical analysis and simulations, optimal pitches that result in local minima of the ripple amplitude with respect to the jaw width and off-axis distance were identified and verified. The effectiveness of optimization in reducing the thread effect were also studied. RESULTS Analysis and simulation based on the square-shaped jaw profiles well characterize the thread effect. Simulations based on the real jaw profiles show reduced ripples and very good agreement of optimal pitches compared with those based on the square profiles. The optimal pitches were found to have little jaw width dependence, except for the real jaw profile of small width (1.05 cm). The optimal pitches for the real jaw profile of width 1.05 cm are unidentifiable except for the largest ones, due to the relative smoothness of the jaw profile. With optimized intensity modulation, the thread effect can be largely suppressed. For real jaw profiles, the optimal pitches with or without dose optimization do not change much. The numbers 0.86∕n found by Kissick et al. well approximate the optimal pitches for off-axis distance of 5 cm. However, optimal pitches are not universal for different off-axis distances: they decrease as the off-axis distance increases. CONCLUSIONS The thread effect can be well explained by the proposed model. Optimization can largely reduce the thread effect. However, an optimal pitch reduces the ripple much easier especially when optimization is limited by many constraints. The optimal pitches predicted by the proposed model could be used as a reference for pitch selection regardless the tumor is at large or small off-axis distance.

Adaptive fractionation therapy: II. Biological effective dose

Radiation therapy is fractionized to differentiate the cell killing between the tumor and organ at risk (OAR). Conventionally, fractionation is done by dividing the total dose into equal fraction sizes. However, as the relative positions (configurations) between OAR and the tumor vary from fractions to fractions, intuitively, we want to use a larger fraction size when OAR and the tumor are far apart and a smaller fraction size when OAR and the tumor are close to each other. Adaptive fractionation accounts for variations of configurations between OAR and the tumor. In part I of this series, the adaptation minimizes the OAR (physical) dose and maintains the total tumor (physical) dose. In this work, instead, the adaptation is based on the biological effective dose (BED). Unlike the linear programming approach in part I, we build a fraction size lookup table using mathematical induction. The lookup table essentially describes the fraction size as a function of the remaining tumor BED, the OAR/tumor dose ratio and the remaining number of fractions. The lookup table is calculated by maximizing the expected survival of OAR and preserving the tumor cell kill. Immediately before the treatment of each fraction, the OAR-tumor configuration and thus the dose ratio can be obtained from the daily setup image, and then the fraction size can be determined by the lookup table. Extensive simulations demonstrate the effectiveness of our method compared with the conventional fractionation method.

Patient-specific dosimetric endpoints based treatment plan quality control in radiotherapy.

In intensity modulated radiotherapy (IMRT), the optimal plan for each patient is specific due to unique patient anatomy. To achieve such a plan, patient-specific dosimetric goals reflecting each patients unique anatomy should be defined and adopted in the treatment planning procedure for plan quality control. This study is to develop such a personalized treatment plan quality control tool by predicting patient-specific dosimetric endpoints (DEs). The incorporation of patient specific DEs is realized by a multi-OAR geometry-dosimetry model, capable of predicting optimal DEs based on the individual patients geometry. The overall quality of a treatment plan is then judged with a numerical treatment plan quality indicator and characterized as optimal or suboptimal. Taking advantage of clinically available prostate volumetric modulated arc therapy (VMAT) treatment plans, we built and evaluated our proposed plan quality control tool. Using our developed tool, six of twenty evaluated plans were identified as sub-optimal plans. After plan re-optimization, these suboptimal plans achieved better OAR dose sparing without sacrificing the PTV coverage, and the dosimetric endpoints of the re-optimized plans agreed well with the model predicted values, which validate the predictability of the proposed tool. In conclusion, the developed tool is able to accurately predict optimally achievable DEs of multiple OARs, identify suboptimal plans, and guide plan optimization. It is a useful tool for achieving patient-specific treatment plan quality control.

A slit method to determine the focal spot size and shape of TomoTherapy system

PURPOSE To obtain accurate x-ray source profile measurements using a slit-collimator, the slit-collimator should have a narrow width, large height, and be positioned near the source. However, these conditions may not always be met. In this paper, the authors provide a detailed analysis of the slit measurement geometry and the relationship between the slit parameters and the measured x-ray source profile. The slit model allows the use of a shorter and more easily available slit-collimator, while accurate source profile measurements can still be obtained. METHODS Measurements were performed with a variety of slit widths and/or slit to source distances. The relationship derived between the slit parameters and the measured profile was used to determine the true focal spot profile through a least square fit of the profile data. The model was verified by comparing the predicted profiles at a variety of slit-collimator parameters with the measured results on the TomoTherapy Hi-Art system. RESULTS Both the treatment beam and the imaging beam were measured. For treatment mode, it was found that a source consisting of one Gaussian with a 0.75 mm full-width-half-maximum (FWHM) and 72% peak amplitude and a second Gaussian with a 2.27 mm FWHM and 18% peak amplitude matched measurement profiles. The overall source profile has a FWHM of 0.93 mm, but with a higher amplitude in the tail region than a single Gaussian. For imaging mode, the source consists of one Gaussian with a 0.68 mm FWHM and 82% peak amplitude and a second Gaussian with a 1.83 mm FWHM and 18% peak amplitude. The overall source profile has a FWHM of 0.77 mm. CONCLUSIONS Our study of the focal spot measurement using slit-collimators showed that accurate source profile measurements can be achieved through fitting of measurement results at different slit widths and source-to-slit distances (SSD). Quantitative measurements of the TomoTherapy linac focal spot showed that the source distribution could be better described with a model consisting of two Gaussian components rather than a single Gaussian model as assumed in previous studies.

Quantitative characterization of tomotherapy MVCT dosimetry.

Megavoltage computed tomography (MVCT) is used as image guidance for patient setup in almost every tomotherapy treatment. Frequent use of ionizing radiation for image guidance has raised concern of imaging dose. The purpose of this work is to quantify and characterize tomotherapy MVCT dosimetry. Our dose calculation was based on a commissioned dose engine, and the calculation result was compared with film measurement. We studied dose profiles, center dose, maximal dose, surface dose, and mean dose on homogeneous cylindrical water phantoms of various diameters for various scanning parameters, including 3 different jaw openings (of nominal value J4, J1, and J0.1) and couch speeds (fine, normal, and coarse). The comparison between calculation and film measurement showed good agreement. In particular, the thread pattern on the film of the helical delivery matched very well with calculation. For the J1 jaw and coarse imaging mode, the maximum difference between calculation and measurement was about 6% of the center dose. Calculation on various sizes of synthesized phantoms showed that the center dose decreases almost linearly as the phantom diameter increases, and that the fine mode (couch speed of 4mm/rotation) received twice the dose of the normal mode (couch speed of 8mm/rotation) and 3 times that of the coarse mode (couch speed of 12mm/rotation) as expected. The maximal dose ranged from 100% to ∼200% of the center dose, with increasing ratios for larger phantoms, smaller jaws, and faster couch speed. For all jaw settings and couch speeds, the mean dose and average surface dose vary from 95% to 125% of the center dose with increasing ratios for larger phantoms. We present a quantitative dosimetric characterization of the tomotherapy MVCT in terms of scanning parameters, phantom size, center dose, maximal dose, surface dose, and mean dose. The results can provide an overall picture of dose distribution and a reference data set that enables estimation of CT dose index for the tomotherapy MVCT.