US6792073B2  Method and apparatus for radiotherapy treatment planning  Google Patents
Method and apparatus for radiotherapy treatment planning Download PDFInfo
 Publication number
 US6792073B2 US6792073B2 US09850533 US85053301A US6792073B2 US 6792073 B2 US6792073 B2 US 6792073B2 US 09850533 US09850533 US 09850533 US 85053301 A US85053301 A US 85053301A US 6792073 B2 US6792073 B2 US 6792073B2
 Authority
 US
 Grant status
 Grant
 Patent type
 Prior art keywords
 dose
 distribution
 carlo
 monte
 denoising
 Prior art date
 Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
 Expired  Fee Related, expires
Links
Images
Classifications

 A—HUMAN NECESSITIES
 A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
 A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
 A61N5/00—Radiation therapy
 A61N5/10—Xray therapy; Gammaray therapy; Particleirradiation therapy
 A61N5/103—Treatment planning systems
 A61N5/1031—Treatment planning systems using a specific method of dose optimization

 A—HUMAN NECESSITIES
 A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
 A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
 A61N5/00—Radiation therapy
 A61N5/10—Xray therapy; Gammaray therapy; Particleirradiation therapy
 A61N5/103—Treatment planning systems
 A61N5/1031—Treatment planning systems using a specific method of dose optimization
 A61N2005/1034—Monte Carlo type methods; particle tracking
Abstract
Description
The present application claims priority from provisional application No. 60/201,959, filed May 5, 2000, entitled “Method of Radiotherapy Treatment Planning”, which is incorporated in its entirety herein by reference.
This invention relates generally to accelerating Monte Carlo calculations, and more particularly to determining radiation dose distributions in a patient for radiotherapy treatment planning by accelerating Monte Carlo dose distribution calculations through denoising of raw Monte Carlo dose distributions.
Medical equipment for radiation therapy treats tumorous (malignant) tissue with high energy radiation. The amount of radiation and its placement must be accurately controlled to ensure that the tumor receives sufficient radiation to be destroyed and that damage to the surrounding and adjacent nontumorous tissue is minimized.
Internal source radiation therapy (referred to as brachytherapy) places capsules of radioactive material inside the patient in proximity to the tumorous tissue. Dose and placement are accurately controlled by the physical positioning of the isotope. However, internal source radiation therapy has disadvantages similar to those present with surgically invasive procedures, which include patient discomfort and risk of infection.
External source radiation therapy uses high energy radiation, collimated to direct a beam into the patient to the tumor site. Although the size and strength of the radiation beam from the external source may be accurately controlled outside of the patient, the dose received by a given volume within the patient may vary because of radiation scattering and absorption by intervening tissue. For these reasons, a determination of the proper dose and placement of the dose requires an estimation of the effects of treated tissue and the tissue surrounding the treated area in scattering or attenuation of the radiation beam.
The Monte Carlo method is the most accurate method for predicting dose distributions. The path of the particles (electrons, photons, protons, or neutrons) through the patient are simulated by Monte Carlo software. Particle interactions (scattering, attenuation, and energy deposition) are simulated one at a time. The Monte Carlo method traces paths of several million particles through a patient model, the patient model accurately reflecting the three dimensional variations of electron density within the volume of the patient under study. For large numbers of source radiation particles (typically above 10^{7}), the Monte Carlo method produces an accurate representation of the dose distribution. For these reasons, the Monte Carlo method is clinically preferred for the calculation of radiation dose in electron beam radiotherapy.
Unfortunately, the Monte Carlo method is extremely time consuming, taking upwards of an hour to compute a single dose distribution using current computers. One crucial and unique feature to Monte Carlo results is that the Monte Carlo method has no welldefined preset “finish” time. Instead, dose distributions can be produced at any time, even after simulating only one source particle. However, the resulting dose distributions are “noisy” (distorted), in inverse proportion to the number of source particles simulated. The Monte Carlo calculations are terminated when the noise level falls below a level deemed acceptable by the user.
Often, radiotherapy treatment planners wish to compare many dose distributions before selecting a final distribution for treatment. Hence, the need exists for a method of dose modeling which is as accurate as the Monte Carlo method but which has greater computational efficiency than the Monte Carlo method alone. This need especially exists as Monte Carlo programs are being developed for most commercial clinical systems.
More generally, the problem of noise is inherent in all Monte Carlo calculations that produce values on a grid. The problem of noise exists whether the grid elements are regularly spaced or irregularly spaced. For these reasons, there also exists a need to reduce noise and provide an accurate and computationally efficient estimate for Monte Carlo calculations in general. The present invention also satisfies this need, providing a method hereafter referred to as Monte Carlo denoising.
The present invention provides an estimate of an actual radiation dose distribution by reducing noise (denoising) from raw results of Monte Carlo generated dose distributions. The resulting denoised dose distribution more quickly and efficiently converges to the required level of accuracy for radiotherapy treatment planning than Monte Carlo calculations alone, and at a reduced computational cost. Optimal denoising reduces Monte Carlo run times by a factor of at least 3 to 6.
The present invention provides a method of accounting for details of the radiation source geometry and materials, and local changes in density of the patient at different points in the irradiated volume. This is done by denoising Monte Carlo dose distributions without significantly increasing the running time or introducing distortions to the underlying dose distributions.
The insight underlying the present invention is that a true (noiseless) radiation dose distribution is smoother (more spatially coherent) than a raw Monte Carlo result. In the present invention, it is the dose distribution itself that is accurately smoothed to reduce noise while not distorting the true underlying (signal) dose distribution.
It is our novel insight into the details of radiation transport physics which indicates that denoising is feasible and desirable (Deasy, June 2000, Denoising of electron beam Monte Carlo dose distributions using digital filtering techniques, Physics in Medicine and Biology, 45:17651779), which publication is incorporated in its entirety herein by reference. In simple terms, diffusive radiation transport implies that the actual dose distribution is smooth, whereas the Monte Carlo result without denoising contains statistical fluctuations which are more rough than the expected underlying dose distribution. Therefore, appropriate denoising techniques can be used to reduce the higher frequency noise of Monte Carlo results while not distorting the true underlying (signal) dose distribution.
In one aspect of the present invention, a method and an apparatus for accomplishing the method is provided where a raw Monte Carlo data distribution D0 (D0 everywhere refers to the raw Monte Carlo result) is first obtained and then denoised to produce a denoised distribution D1 (D1 everywhere refers to the denoised Monte Carlo result). This aspect of the present invention could be used to simply accelerate the computation of any Monte Carlo result or the present invention could be more specifically used to determine a dose distribution in a patient for radiation therapy treatment planning.
In another aspect of the present invention, denoising of the raw Monte Carlo dose distribution D0 could employ digital filtering, wavelet denoising, kernel smoothing or nonparametric regression smoothing. The digital filtering techniques could employ Binomial/Gaussian filters or localleast squares filters.
In another aspect of the present invention, denoising the raw Monte Carlo dose distribution D0 could first involve transforming the Monte Carlo dose distribution D0 by computing the square root of each data element of the Monte Carlo dose distribution D0. Digital filtering techniques could then be applied to the transformed data elements to obtain a filtered result. Then, the elements of the filtered result are squared to obtain a final best estimate of the dose distribution. The digital filtering techniques could employ Binomial/Gaussian filters or localleast squares filters. If employed, the Binomial/Gaussian filters could use a higher degree polynomial for less aggressive denoising. If employed, the localleast squares filters could use a smaller neighborhood for less aggressive denoising.
In a further aspect of the present invention, denoising the raw Monte Carlo dose distribution D0 employs either a spatially adaptive iterative filtering (SAIF) algorithm, a wavelet shrinkage threshold denoising algorithm, a spatially adaptive wavelet denoising (SAWD) algorithm, or a batchaveraged wavelet denoising (BAWD) algorithm.
In a still further aspect of the present invention, a method for determining a radiation dose pattern directed at a volume of interest in a patient includes providing to a computer a matrix T (describing an electron density for elements of the volume of interest in the patient), a subroutine (to produce a list of source particles describing a radiation source) and a Monte Carlo program (to simulate effects of the source particles). A dose accumulation matrix A is initialized to zeros. The computer then runs the Monte Carlo program to produce a dose deposition pattern for an incident source particle. The dose accumulation matrix A is then updated by adding the dose deposition pattern produced for the incident source particle. Running the Monte Carlo program and updating the dose accumulation matrix A is repeated until a predetermined stopping criteria is reached. The Monte Carlo dose distribution D0 is then set equal to the updated dose accumulation matrix A. Finally, a denoised estimate is developed for each element of the Monte Carlo data distribution D0 to produce a denoised dose distribution D1 representative of the radiation dose pattern. This aspect of the present invention could be varied by selecting an incident source particle to reduce uncertainty in a predetermined location of a dose deposition pattern before running the Monte Carlo program. Selecting the incident source particle then becomes part of the repeated process (along with running the Monte Carlo program and updating the dose accumulation matrix A) until the predetermined stopping criteria is reached.
In a yet further aspect of the present invention, developing the denoised estimate includes selecting an element from the Monte Carlo data distribution D0, creating a submatrix using values from the Monte Carlo data distribution D0 within a predetermined range of the selected element and multiplying the submatrix by a predetermined matrix designed to effect a leastsquares fit then extracting the fit to the center element of the submatrix as the leastsquares fit, and repeating these acts for each element of the Monte Carlo data distribution D0 to complete the development of the denoised estimate.
In another aspect of the present invention, developing the denoised estimate uses a variance stabilizing transformation (defined as a transformation which results in the statistical having more similar characteristics in every part of the data distribution). The variance stabilizing transformation could include computing the square root of each element of the Monte Carlo data distribution D0 to create a corresponding matrix D2. The corresponding matrix (D2) is then denoised to produce a matrix M. Every element of the matrix M is then squared to produce the denoised estimate.
In still another aspect of the present invention, developing the denoised estimate uses a statisticallybased smoothing method. The statisticallybased smoothing method could be the Loess method.
In yet another aspect of the present invention, developing the denoised estimate uses a curve fitting method employing Hermite polynomials for the leastsquares fit. Developing the denoised estimate could also use a basis function selected from the group consisting of spline fits, spline fits with knots and weighted local least squares methods. Or, developing the denoised estimate could use a Fourierbased smoothing method. The Fourierbased smoothing method could employ a lowpass filter or a Wiener filter.
In a further aspect of the present invention, the Monte Carlo data distribution D0 is divided into parallel planes and developing the denoised estimate occurs independently for each plane.
In a still further aspect of the present invention, the Monte Carlo data distribution D0 is regridded onto a noncartesian grid and developing the denoised estimate occurs on the noncartesian grid with the resulting denoised dose distribution D1 being transformed back to a cartesian grid. The noncartesian grid could be a diverging fan beam grid.
In a yet further aspect of the present invention, developing the denoised estimate uses a variable computational technique that varies to address selected portions of an image.
In another aspect of the present invention, developing the denoised estimate uses a computational technique that reduces random noise and accurately estimates underlying signal.
It is therefore one object of the invention to provide a method of dose calculation that is more computationally efficient at reaching the same level of noise than using present Monte Carlo methods alone to simulate radiation dose distributions.
It is another object of the invention to provide a method of dose calculation which reaches a lower level of dosage pattern noise than Monte Carlo methods for the same computation time.
It is a further object of the invention to accelerate any type of Monte Carlo calculation (for any purpose) to accurately represent a true underlying signal by generating and denoising a raw Monte Carlo calculation, the present invention reducing computation time relative to using present Monte Carlo methods alone.
For the purpose of illustrating the invention, there is shown in the drawings a form which is presently preferred; it being understood, however, that this invention is not limited to the precise arrangements and instrumentalities shown.
FIG. 1 illustrates a digital stepheterogeneity phantom, where electrons are transported through air before reaching the water surface, that is used as input for Monte Carlo calculations to demonstrate the denoising techniques of the present invention;
FIG. 2 illustrates profile plots demonstrating digital filtering techniques to a uniform water phantom in accordance with the present invention;
FIG. 3 illustrates profile plots demonstrating digital filtering techniques to the stepheterogeneity phantom in accordance with the present invention;
FIG. 4 illustrates contour plots of dose distributions resulting from the uniform water phantom and the heterogeneity phantom in accordance with the present invention;
FIG. 5 illustrates application of spatially adaptive iterative filtering (SAIF) denoising algorithm to the heterogeneity phantom in accordance with the present invention;
FIG. 6 illustrates denoising of the heterogeneity phantom using wavelet denoising algorithm with 9,7 biorthogonal filters in accordance with the present invention;
FIG. 7 illustrates denoising of the heterogeneity phantom using wavelet shrinkage denoising with thirdorder ‘Coiflet’ wavelets and soft thresholding in accordance with the present invention; and
FIG. 8 illustrates wavelet denoising of threedimensional photon dose distributions, in accordance with the present invention.
The present invention determines dosage patterns in radiotherapy where internal or external radiation sources are used to direct radiation at a volume of interest in a patient. The resulting dosage pattern within the patient is first deduced a Monte Carlo simulations tracing the path and energy deposition events of radiation particles emitted by the source. The energy deposition events are stored in memory locations associated with defined cuboid volumes of interest within the patient (voxels).
The method employs an electronic computer into which a matrix T is entered, the matrix T describing patient tissue density at the elements within the volume of interest of the patient. Radiation source particles are introduced to the matrix T, the radiation source particles having characteristics statistically consistent with the source as determined from other sources of information well known to practitioners in the art.
The dose pattern is initially estimated using Monte Carlo methods also well known to practitioners of the art. The dose to each voxel in the dose distribution is then estimated more accurately by reducing the noise in the Monte Carlo dose distribution, thereby producing a more accurate estimate of the actual smooth radiation dose pattern.
The Monte Carlo Method and Noise Characteristics
Variance at each voxel (memory location associated with a defined volume within a patient) varies approximately linearly with the corresponding absolute Monte Carlo dose. A Monte Carlo dose distribution can be considered a dose “image” measured by perfect “voxel detectors” which sum the energy deposition events. It is as though a random Gaussian variable is added to each squareroot dose value, where the width of the random variable's probability density function is proportional to the dose value to which it is added. The squareroot transform (i.e. taking the squareroot of all dose values) is the appropriate variance stabilizing transform for such statistics (Poisson statistics). Hence the model:
Using the model, the new array (whose values equal the square root of the dose array) consists of the signal or image to be recovered (True Dose Array^{1/2}), plus an additive random variable noise term of mean value zero and slowly varying—nearly constant—standard deviation.
The data model of Equation 1 of this section is known as a stationary (not a function of position) Gaussian noise model. This is the most basic data models in general technical use and it arises in many different technical areas, such as statistics, digital filtering, and digital image processing.
Denoising of Monte Carlo Electron Beam Dose Distributions via Digital Filtering
In various aspects of the present invention, four different filters (FIG. 2) are applied to four ‘benchmark’ style MONTE CARLO results to examine the effect of filtering on resulting dose distributions. In one embodiment the invention uses twodimensional filtering and image adaptive methods, although it will be apparent to those skilled in the art that the present invention is equally applicable to threedimensional filtering and image adaptive methods. The filtering techniques are applied to a variance stabilized (squareroot) dose image which is transformed back to dosespace (i.e. squared) after filtering.
Denoising is aimed at computing a new image h, having elements h_{i,j}, such that
where the a_{m,n }values are the filter coefficients, the g values are the original MC dose values, and the local convolution (smoothing) is a function of dose values within a distance ±k. The filters applied are listed in Table 0. The filter classes are of two types: Binomial/Gaussian filters (lowpass filters) and localleast squares (SavitzkyGolay) filters.
TABLE 0  
Listing of digital filters. The ‘order’ of each Gaussian filter  
(GAUSS) is the number of applications of a nearestneighbor  
binomial averaging that is used to derive that filter. The ‘order’  
of each local leastsquares (LLS) filter is the highest polynomial  
power from which that filter's coefficients are derived  
Filter name  Description 
GAUSS_W5_O4  Gaussian (binomial) filter of width 5 and order 4 
GAUSS_W3_O2  Gaussian (binomial) filter of width 3 and order 2 
LLS_W11_O5  localleastsquares filter of width 11 and order 5 
LLS_W5_O3  localleastsquares filter of width 5 and order 3 
Binomial/Gaussian filters are constructed in 1D by repeatedly creating a new vector in which each bin is ½ the sum of nearest neighbors in the original vector of bin values. The process begins with a single nonzero bin equal to one. After each step, the new bins are centered at the borders of the old bins. The order of the filtering process is equal to the number of times this process is repeated.
In 2D, the operation is performed separately on each row of the filter matrix, and then on each column; this cycle is repeated a number of times equal to the filter order. The filter is applied by convolution with the original image. Binomial/Gaussian filters result in diffusionlike behavior and a Gaussianlike Fourier transfer functions. Hereafter, these filters are simply referred to as Gaussian filters.
Another interesting class of filters arises from the application of local leastsquares principles. The neighborhood around each data point is fit using a polynomial fitting function, under the assumption that each data point is characterized by the same variance, according to the equation:
where unit distance is fixed as the voxel width. The order of the filter is the order of the polynomial. This equation is used to find an estimate of the center voxel in the neighborhood, which then replaces that voxel's value in the new denoised image. This process is repeated such that every voxel becomes in turn the center voxel in the application of Equation (3).
This type of filter is called a Local LeastSquares filter (LLS) and is also known as a SavitzkyGolay filter in applications to chemical spectral analysis. Previous applications of LLS filters to 2D data are rare. We developed Mathematica code to generate LLS coefficients for an arbitrary filter width, fitted polynomial order, and data dimension (1, 2, or 3D).
FIG. 2 shows filtered results for the uniform water phantom at a depth of 1.0 cm. Six profile plots, through the center of the dataset (column 60) are compared with the reference high density result, which is always shown as the thick line in all figures. Voxel midpoints match step midpoints in all figures. The upper pair of plots compare the high and low density results; the other four plots compare the high density result with the denoised low density result foil the indicated filter. The voxel widths are 0.2 cm. In the plateau, all the filters improve substantially on the low density result. However, there is some blurring of the sharp edges. FIG. 3 shows filtered results for the stepheterogeneity phantom of FIG. 1, are through data column 60 (one of the two centermost columns in the 2D dose matrix), are always compared with the reference high density result (thick lines), and are at a depth of 2.55 cm below the top of the step heterogeneity. The low density dose standard deviation is 2.9%. This is an extreme test of possible heterogeneity effects. The filtering process clearly reduces noise, although some blurring of the dipandrise is evident for all filters except LLS_W5_O3, where it cannot be visually detected. The blurring is worst for Gauss_W5_O4.
FIG. 4 shows contour plots of the dose distributions resulting from the uniform water phantom (left column) and the heterogeneity phantom (right column). The left column is from homogeneous water phantom results at a depth of 2.5 cm; the right column is from the heterogeneous phantom at a depth of 2.55 cm below the top of the step heterogeneity. These plots emphsize the difficulty in obtaining a truly ‘clean’ contour plot relying on ‘raw’ Monte Carlo results. For both phantoms, the denoised distribution (using LLS_W5_O3) is clearly cleaner visually than the low density Monte Carlo simulation. The bottom row shows the result of filtering the low density Monte Carlo data.
Denoising using digital filtering has the following advantages:
1. Improved visualization of isodose lines (FIG. 3). The ability to visually inspect the dose distribution is significantly improved using denoising filtering.
2. Suppression of spurious dose excursions (noise), as seen in the profile plots (FIGS. 23).
There is a tradeoff between local denoising aggressiveness and systematic distortion (bias) due to blurring, as might be expected based on frequency considerations. Nonetheless, for electron beams, the filter LLS_W5_O3 is a candidate for routine filtering of clinically produced Monte Carlo dose distributions.
Spatially Adaptive Iterative Filtering (SAIF)
In another aspect of the present invention, a oneparameter denoising algorithm is employed and is directed to locally filter as aggressively as possible without distorting the underlying dose distribution. This spatially adaptive iterative filtering (SAIF) algorithm uses an array of denoised dose images by applying an array of digital filters to the raw Monte Carlo image. The best denoised image is constructed by selecting, for each voxel, the most aggressive denoising result which is statistically consistent with the original dose image. The statistical test is a chisquare comparison, over a local region, between the digital filter result and a ‘lightly’ denoised result which is known to be not too aggressive.
The method of spatially adaptive iterative filtering (SAIF):
Step 1. Begin with a safely denoised image (no significant bias), using the LLS_W5_O3 filter. The result is denoted F_{0}.
Step 2. Set the current best dose image (CBD) equal to F_{0}.
Step 3. Loop over increasingly aggressive filters:
Step 3.1. Denoise with the i'th filter, yielding F_{1}.
Step 3.2. Form a matrix of ‘local chisquare values,’ where each value is the chisquare comparison between the F_{1 }doses and the F_{0 }doses over a small region centered about that voxel.
Step 3.3. For each position in the chisquare image:
Step 3.3.1. If the local chisquare is below a preset threshold, replace that position in CBD with dose at that position in F_{1}.
FIG. 5 illustrates an example using SAIF. The local chisquare test determines if the new aggressive filter ‘fits’ the image. In the FIG. 5 embodiment, an application of a (SAIF) algorithm to the electron beam heterogeneity phantom is shown. In the middle view, profile plots across the heterogeneity, at column 55, are shown. Successive results with progressively mote aggressive filters are given in the 3rd, 4th and 5th subplots. The fifth subplot is the final result. In this embodiment, dose estimates from the two more aagressive filters are only accepted if a chisquare value comparing each result and the safe (mild filter) result, over a local 3×3 region, is less than 0.0001 (the dose images are normalized to a maximum value of 1). This threshold is chosen empirically on the basis of the chisquare histogram shown in the second panel. The histogram shows how the local chisquare value is sensitive to bias which results in the two extra peaks shown. The top panel shows the difference between the resulting dose values and the “true” dose values (144 million source electrons in this case). The lack of any clear pattern in the difference indicates that the underlying smooth dose distribution is wellapproximated by the SAIF algorithm.
Wavelet Denoising
In a further aspect of the present invention, wavelet shrinkage threshold denoising, or simply wavelet denoising is employed to denoise the raw Monte Carlo dose distribution. Under a wavelet transform, the data set is expressed as a linear combination or expansion in discrete wavelet coefficients. Typically, small coefficient values in the expansion represent noise, whereas larger coefficients represent spatially coherent structures. After transforming the data, the coefficients below some threshold are set to zero (in ‘hard thresholding’). Then, reverse wavelet transformation is performed which converts the truncated array of wavelet coefficients back to dose array values.
In one embodiment, a Wavelet Shrinkage Denoising Algorithm with 9,7biorthogonal filters is used, the process summarized as follows:
Step 1. Apply the squareroot variance stabilizing transformation.
Step 2. Spincycling: Create 8 extra matrices by shifting the original by one pixel value in all the nearest neighbor directions. (This eliminates grid orientation artifacts.)
Step 3. On the original plus the 8 shifted arrays do the following:
Step 3a. Transform: perform the forward wavelet transform using the 9,7biorthogonal filters.
Step 3b. Hard thresholding: set all wavelet coefficients smaller than a fixed fraction (denoted epsilon) of the maximum wavelet coefficient equal to zero (those coefficients are predominantly noise).
Step 3c. Inverse transform: perform the 9,7 inverse wavelet transform.
Step 4. Backshift the 8 shifted arrays.
Step 5. Average the backshifted arrays and the unshifted array.
Step 6. Inverse variance stabilization: square the resulting array values (step). This is the denoised dose distribution.
FIG. 6 illustrates denoising of the heterogeneity dose distribution of FIG. 1 using the wavelet denoising algorithm with 9,7biorthogonal filters. Wavelet coefficient values less than 0.1% of the maximum wavelet coefficient value were set to zero. A computationally efficient liftingscheme was implemented in ANSI C. The profile is through the center of the dose distribution (column 60). The thick line is the result of wavelet shrinkage denoising applied to a lowdensity result (thin line), which had a dose standard deviation of 2.9% before denoising. The computer time required for the denoising of this 120×120 matrix is a few milliseconds. At 0.01 threshold level, there is no perceptible bias even though noise is significantly reduced over much of the profile plot. One advantage of wavelet methods over traditional smoothing techniques is that they are locally adaptive in an automatic fashion.
FIG. 7 illustrates an example of wavelet denoising using a 3rd order ‘Coiflet’ wavelet and ‘soft thresholding’. Soft thresholding, like hard thresholding, sets all wavelet coefficients below a threshold to zero, but also ‘shrinks’ coefficients above the threshold toward zero by an amount equal in magnitude to the threshold. In the FIG. 7 embodiment, profile plots are for the stepheterogeneity phantom at a depth of 2.55 cm. Plots are through column 60 (one of the two centermost columns in the 2D dose matrix), and are compared with the reference high density result (thick lines). The different subplots represent thresholds which have been scaled by factors of 10. The low density dose standard deviation is 2.9%. The lower left figure utilizes the optimal (or nearly optimal) threshold value.
FIG. 8 illustrates wavelet denoising for 3D photon dose distributions. The Pencil beam (PB) is denoised by denoising 2D slices of the dose distribution in the xz plane (i.e., a transverse“beam direction” plane) using the 9,7biorthogonal code. Little bias is introduced and run time is effectively reduced by at least a factor of four. In the FIG. 8 embodiment, ICCR phantom geometry was used which consists of slabs of water (03 cm depth), aluminum (35 cm), homogeneous lunglike media (512 cm), and water (1232 cm). The NRCC 6 MV spectrum specified by the ICCR test is used in the ITS code. The beam is collimated to 2 cm×2 cm by an upstream lead collimator. The top shows a Grayscale central slice, the next down is a highdensity central axis depthdose with 1.2 billion source particles, the next shows a central axis depthdose with 280 million source particles denoised with a 0.25% threshold, and the bottom illustrates a relative difference between the two. Although there may be some systematic differences, the high density result is noisier than the low density result, indicating that speed gain is greater than a factor of four.
A thorough description of wavelet denoising using 9,7biorthogonal filters follows: The simulated dose distribution d consists of two parts, d=s+n, where s is the smooth function obtained in the limit by running the Monte Carlo simulation forever; and n is the rough function giving error due to the short Monte Carlo simulation. The data d is sampled in space, giving an array of values that are then linearly transformed into discrete wavelet coefficients, described as W(d)=W(s+n)=W(s)+W(n).
A key property of the wavelet transform is that a sampled rough function n gives values W(n) which are more nearly equal in amplitude than those of a sampled smooth function s. If n is relatively small compared to s, then any sufficiently small wavelet coefficient is more likely be part of n than s. By picking a threshold δ>0 and setting W(d)_{i}=0 if W(d)_{i}<δ, attenuate the n, or noise component of d. The array reconstructed from the surviving coefficients is therefore a closer approximation to s than was d=s+n.
For the key property to hold; it is necessary that the wavelet transform use a wavelet that is about as smooth as s, and thus smoother than n. The 9,7biorthogonal wavelet is chosen, as used by the draft JPEG 2000 image compression algorithm, since images are about as smooth as true underlying dose distributions.
Approximations of Dose Distributions
Monte Carlo radiotherapy simulations produce sampled dose distributions on a finite grid. A simulated twodimensional dose distribution on an M×N grid may be represented by a nonnegative function f=f(m,n), where 0≦m<M and 0≦n<N are integer indices. Such functions may be considered approximations to distributions φ(x,y) defined on a continuum 0≦x,y<1 as a superposition of basic distributions φ_{mn}=(x,y) concentrated near
given by the formula
For example, characteristic functions of small regions could be taken:
A superposition of basic densities φ^{0} _{mn }is constant on
rectangles of the unit square. That is not a smooth or even continuous function.
Normal probability densities
of suitable mean μ and variance σ, give smooth basic functions:
A superposition of densities φ^{∞} _{mn }is not confined to the unit square, or any bounded region, but it has an attractive physical interpretation if the assumption is made that errors of dose measurement are normally distributed.
The represented φ will be as smooth as its basic densities, regardless of f, so a selection between φ^{0} _{mn }and φ^{∞} _{mn }is made. Basic densities with two continuous derivatives that are “almost” confined to the unit square are selected.
Efficient Coding of Fluctuations as Wavelets
Regions of nearly constant dosage might span many grid points, with relatively small fluctuations on various scales. The amplitude and scale of each fluctuation can be used to decide what is noise and what is signal.
Independent wave function fluctuations of different sizes and positions are modeled. By approximating the original dose distribution as a superposition of particular basic densities, the fast discrete wavelet transform algorithm with 9,7biorthogonal symmetric wavelets can be used for computing the scale and amplitude of fluctuations.
The symmetry of boundary wavelets breaks when the dimensions are not powers of 2, but both the wavelets and scaling functions remain smooth functions when reflected at the right boundary and periodized.
For efficiency, the lifting scheme is used. To reduce bias arising from the arbitrary choice of grid origin, the results of denoising the original image, and the 8 periodic shifts by one grid point: left, right, up, down, and four diagonals, are averaged. This is necessary when using the 9,7biorthogonal symmetric wavelets, but is not necessary if other wavelet families are used as the basis of the denoising.
Discrete Wavelet Transforms
The wavelet transform is identical to that of JPEG 2000. For twodimensional and threedimensional arrays, the transform is applied separately along each dimension, so only the onedimensional version is described. Let x=x(k), 0≦k<K be the array of K values to be transformed. Generate a newarray Fx=Fx(k), 0≦k<K, by the following rule:
Here h=h(n) and g=g(n) are the filter sequences defining the 9,7biorthogonal wavelets. Table 1 gives their approximate values. The index ranges [−4, 4] and [−4, 2] have 9 and 7 indices, respectively. Before summing, extend the original signal x to the longer signal X by wholesample symmetric reflection, as
TABLE 1  
9,7biorthogonal analysis filters h, g and their inverse, or  
synthesis filters h′, g′.  
k  h(k) = (−1)^{k}g′(k)  k  g(k) = −(−1)^{k}h′(k) 
−4, 4  0.03782845550699  −4, 2  0.064538882628938 
−3, 3  −0.02384946501938  −3, 1  −0.040689417609558 
−2, 2  −0.11062440441842  −2, 0  −0.418092273222212 
−1, 1  0.37740285561265  −1  0.788485616405664 
0  0.85269867900940  
defined by the formula below
(Note: the above table and equation could use a higher resolution capture, but this is not absolutely necessary)
X is then extended periodically by defining X(k±[2K−2])=X(k) for all integer indices k. The result is (2K−2)periodic and symmetric with respect to reflection about indices 0 and K−1, with X(k)=x(k) for the original indices k=0, 1, . . . , K−1. The output sequence Fx(k) is also (2K−2)periodic and defined at all integers k, but because the filters are symmetric, a complete set of output values may be found at just the indices 0, 1, . . . , K−1. This allows computation of Fx in about 8K operations.
The lifting scheme reorganizes the computation of Fx to save operations. With x extended to X as before, F is performed sequentially in place as follows:
1X(k)←X(k)+α[X(k−1)+X(k+1)], for all odd k in the range 0<k<K;
2. X(k)←X(k)+β[X(k−1)+X(k+1)] for all even k in the range 0≦k<K;
3. X(k)←X(k)+γ[X(k−1)+X(k+1)], for all odd k in the range 0<k<K;
4. X(k)←X(k)+δ[X(k−1)+X(k+1)] for all even k in the range 0≦k<K;
5. Fx(k)←X(k)/ζ, for all even k in the range 0≦k<K;
6. Fx(k)←ζX(k), for all even k in the range 0≦k<K;
TABLE 2  
Lifting coefficients for 9,7biorthogonal analysis filters.  
Coefficient  Value  
α  −1.58613434205992  
β  −0.05298011857296  
γ  0.88291107553093  
δ  0.44350685204397  
ζ  1.14960439886024  
The coefficients {α,β,γ,δ,ζ} are given in Table 2. The resulting sequence Fx is the same as the one defined by Equation 1 of this section, but it takes only 5K operations to compute, less intermediate storage, and fewer data exchanges.
To invert F, use {−α,−β,−γ,−δ; 1/ζ} and apply the operations in reverse order:
1. Y(k)←Y(k)/ζ, for all even k in the range 0≦k<K;
2. Y(k)←Y(k)/ζ, for all odd k in the range 0<k<K;
3. Y(k)←Y(k)−δ[Y(k−1)+Y(k+1)], for all even k in the range 0≦k<K;
4. Y(k)←Y(k)−γ[Y(k−1)+Y(k+1)], for all odd k in the range 0<k<K;
5. F^{1}y(k)←Y(k)−β[Y(k−1)+Y(k+1)], for all even k in the range 0≦k<K;
6. F^{1}y(k)←Y(k)−α[Y(k−1)+Y(k+1)], for all odd k in the range 0<k<K.
The analog of Equation 1 of this section for F^{1}y is
where Y is the same extension of y as X is of x. The inverse filters h′, g′ are given in Table 1. The index ranges are K_{0}(n)={k:−4≦2k−n≦2} and K_{1}(n)={k:−4≦2k−n≦4}, for n=0,1, . . . , K−1.
Each evenindexed element Fx(2k) is a weighted average of x over nine grid points near 2k, with the weights being the filter coefficients. It is also the amplitude of the basic distribution shaped, but twice as wide, that best represents the signal near 2k. That original signal differs from this wide base by fluctuations on the scale of two grid points. The amplitudes of those fluctuation wavelets are given by the oddindexed elements Fx(1), Fx(3), . . . , etc.
To compute a onedimensional wavelet transform on a signal x, apply F, retain the oddindexed values, extract the evenindexed values to another array of about half the length, call that new array x, and repeat until x contains just a single number:
1. Let j=1.
2. Compute y(k)=Fx(k), k=0,1, . . . , K−1;
3. Extract D_{j}(k)=y(2k+1) for all indices k with 0<2k+1<K.
4. Replace x(k)=y(2k) for all indices k with 0≦2k<K.
5. Replace K←K/2 if K is even, or K←(K+1)/2 if K is odd. This is the number of elements written to x at step 4.
6. If K>1, then increment j←j+1 and go to step 2. Otherwise, stop.
At the termination of the above algorithm, the discrete wavelet transform of the original signal is stored in the arrays D_{1}, D_{2}, . . . , D_{j}, plus y(0). The values in D_{j }are the amplitudes of fluctuations of the signal at scales of 2^{J }grid points; y(0) contains the average value of the signal. J is the depth of the wavelet transform, the maximum value of j in the algorithm above, and is the least integer such that 2^{j}≧K.
A standard way to index wavelet coefficients so as to preserve both size and position information for the fluctuations is to let W(x)=(y(0); D_{j}, D_{2}, . . . , D_{1}), the concatenation of the fluctuation arrays ill the reverse order of their computation.
Approximation and DeNoising
The root mean square sum is used to measure the size of functions on grids M.
Here M is the number of points in the grid. With this definition, the quantity ∥{circumflex over (f)}−f∥ is the RMS, or root mean square error between estimator {circumflex over (f)} and exact value f, while ∥{circumflex over (f)}−f∥/∥f∥ is the relative RMS error. To say n is relatively small compared to s will mean ∥n∥<<∥s∥.
Define the thresholding function T=T_{∈}[x] for a function x on M by letting a=max{x(m):m∈M} and then setting
Thresholding with small ε is a nonlinear transformation, but is close to the identity in the sense that
for every m∈M.
It remains to choose the threshold. Let ∈ be a fixed positive parameter, let 0∈M be the origin in the grid M, and let 0∈M be some fixed grid point. Designate the wavelet transform of ƒ by W(f)={W(ƒ)_{m}:m∈M}, with W(ƒ)_{0 }being the single average coefficient, previously called y(0). The following algorithm is performed:
1. Translate the function ƒ to ƒ_{o}, so that ƒ_{o}(0)=f(o).
2. Compute the wavelet transform y_{o}←W(ƒ_{o}).
3. Apply ε—thresholding z_{o}(m)←T(y_{o}(m)) for all m≠0.
4. Compute the inverse wavelet transform x_{o}←W^{−1}(z_{o}).
5. Translate the function x_{o }to x, so that x(o)=x_{o}(0).
Let O be the origin of M and its nearest neighbors, then average together all the x's obtained for different choices of o∈O. In the onedimensional case, O={−1, 0, +1}, while in the twodimensional case O is the ninepoint set {−1, 0,+1}x{−1, 0,+1}. Denote by E(f, ε) the average of the x's produced by all the shifts of the original f.
Suppose that d=s+n is an approximation to the dose distribution s. Assume that s is known because the Monte Carlo simulation has been run long enough to approximate s adequately. Obtain several independent d's by selecting short disjoint segments of this long simulation. These are all nonnegative functions on the same grid M, so perform the variancestabilizing squareroot operation and get all estimator for s from
The threshold ∈ that minimizes the RMS error
for each of the short segments d that comprise the long run, is found by a bisection search. The average of these thresholds is an estimator for the best threshold to use in all similarsize simulations.
Threshold Dependence on Run Length
Suppose that a Monte Carlo simulation of length R produces an adequate approximation to the limit dose distribution s. Let ε (r) be the optimal threshold for runs of length r<<R, namely the average of the minimizers of Δ(r) for each short run. Similarly, let be the average of the errors
taken over the R/r runs d of length r that comprise the run that produced s.
Spatially Adaptive Wavelet Denoising Algorithm (SAWD):
In another aspect of the present invention, a spatially adaptive wavelet denoising algorithm is employed, which adapts the local wavelet threshold according to a local measure of goodness of fit. In essence, the most aggressive denoising (which shows no significant bias compared to the safely/unaggressively denoised image) is locally applied.
The Method of spatially adaptive wavelet denoising (SAWD):
Step 1. Begin with a safely denoised image (no significant bias), implemented by wden9_{—}7.c computer code or the like with a threshold of 0.05% for this step. The result is denoted F_{0}.
Step 2. Set the current best dose image (CBD) equal to F_{0}.
Step. 3 Loop over increasing thresholds (for example: [0.1, 0.15, 0.2, 0.25]):
Step 3.1. Denoise with the wavelet denoising wden9_{—}7 software and the i'th threshold, yielding F_{1}.
Step 3.2. Form a matrix of ‘local chisquare values,’ where each value is the chisquare comparison between the F_{1 }doses and the F_{0 }doses over a small region centered about that voxel.
Step 3.3. For each position in the chisquare image:
Step 3.3.1. If the local chisquare is below a preset threshold, replace that position in CBD with the value in F_{i}.
Batchaveraged Wavelet Denoising:
In another aspect of the present invention, a Batchaveraged wavelet algorithm is employed to denoise the raw Monte Carlo dose distribution.
The Method of Batchaveraged wavelet denoising (BAWD):
Step 0. Divide the MC computation into N (say, 10) equal batches. For each batch apply the following wavelet denoising algorithm (implemented in wden9_{—}7.c computer code or the like):
Step 1. Apply the squareroot variance stabilizing transformation.
Step 2. Spincycling: Create 8 extra matrices by shifting the original by one pixel value in all the nearest neighbor directions. (This eliminates grid orientation artifacts.)
Step 3. On the original plus the 8 shifted arrays do the following:
Step 3a. Transform: perform the forward wavelet transform using the 9,7biorthogonal filters.
Step 3b. Hard thresholding: set all wavelet coefficients smaller than a fixed fraction of the maximum wavelet coefficient (denoted epsilon) equal to zero (those coefficients are predominantly noise).
Step 3c. Inverse transform: perform the 9,7 inverse wavelet transform.
Step 4. Backshift the 8 shifted arrays.
Step 5. Average the backshifted arrays and the unshifted array.
Step 6. Inverse variance stabilization: square the resulting array values (step). This is the denoised dose distribution.
Step 7. End loop over batches (go to Step 0 if more batches are required).
Step 8. Return the average of the denoised batch solutions as the final answer.
These and other advantages of the present invention will be apparent to those skilled in the art from the foregoing specification. Accordingly, it will be recognized by those skilled in the art that changes or modifications may be made to the abovedescribed embodiments without departing from the broad inventive concepts of the invention. It should therefore be understood that this invention is not limited to the particular embodiments described herein, but is intended to include all changes and modifications that are within the scope and spirit of the invention as set forth in the claims.
Claims (33)
Priority Applications (2)
Application Number  Priority Date  Filing Date  Title 

US20195900 true  20000505  20000505  
US09850533 US6792073B2 (en)  20000505  20010507  Method and apparatus for radiotherapy treatment planning 
Applications Claiming Priority (1)
Application Number  Priority Date  Filing Date  Title 

US09850533 US6792073B2 (en)  20000505  20010507  Method and apparatus for radiotherapy treatment planning 
Publications (2)
Publication Number  Publication Date 

US20020027971A1 true US20020027971A1 (en)  20020307 
US6792073B2 true US6792073B2 (en)  20040914 
Family
ID=22747995
Family Applications (1)
Application Number  Title  Priority Date  Filing Date 

US09850533 Expired  Fee Related US6792073B2 (en)  20000505  20010507  Method and apparatus for radiotherapy treatment planning 
Country Status (2)
Country  Link 

US (1)  US6792073B2 (en) 
WO (1)  WO2001085252A3 (en) 
Cited By (36)
Publication number  Priority date  Publication date  Assignee  Title 

US20040015073A1 (en) *  20020426  20040122  Michael Schell  Target repositioning error correction filter and method 
US20050201516A1 (en) *  20020306  20050915  Ruchala Kenneth J.  Method for modification of radiotherapy treatment delivery 
US20060085223A1 (en) *  20041006  20060420  Jean Anderson  System and user interface for presenting treatment information 
US20060173725A1 (en) *  20050202  20060803  Michele Abraham  Integrated treatment planning and scheduling system 
WO2007014027A2 (en) *  20050722  20070201  Tomotherapy Incorporated  System and method of remotely analyzing operation of a radiation therapy system 
US20070041495A1 (en) *  20050722  20070222  Olivera Gustavo H  Method of and system for predicting dose delivery 
US20070041499A1 (en) *  20050722  20070222  Weiguo Lu  Method and system for evaluating quality assurance criteria in delivery of a treatment plan 
US20080089480A1 (en) *  20061016  20080417  Oraya Therapeutics, Inc.  Portable orthovoltage radiotherapy 
US20090063110A1 (en) *  20030314  20090305  Transpire,Inc.  Brachytherapy dose computation system and method 
US20090110145A1 (en) *  20071025  20090430  Tomotherapy Incorporated  Method for adapting fractionation of a radiation therapy dose 
US20090234175A1 (en) *  20080311  20090917  Cynthia Maier  System and method for imageguided therapy planning and procedure 
WO2009127747A1 (en)  20080414  20091022  Gmv Aerospace And Defence S.A.  Planning system for intraoperative radiation therapy and method for carrying out said planning 
US7620147B2 (en)  20061213  20091117  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US20100053208A1 (en) *  20080828  20100304  Tomotherapy Incorporated  System and method of contouring a target area 
US20100054413A1 (en) *  20080828  20100304  Tomotherapy Incorporated  System and method of calculating dose uncertainty 
US7792249B2 (en)  20071223  20100907  Oraya Therapeutics, Inc.  Methods and devices for detecting, controlling, and predicting radiation delivery 
US20100228116A1 (en) *  20090303  20100909  Weiguo Lu  System and method of optimizing a heterogeneous radiation dose to be delivered to a patient 
US7801271B2 (en)  20071223  20100921  Oraya Therapeutics, Inc.  Methods and devices for orthovoltage ocular radiotherapy and treatment planning 
US7839972B2 (en)  20050722  20101123  Tomotherapy Incorporated  System and method of evaluating dose delivered by a radiation therapy system 
US20110029329A1 (en) *  20080424  20110203  Koninklijke Philips Electronics N.V.  Dosevolume kernel generation 
US20110122997A1 (en) *  20091030  20110526  Weiguo Lu  Nonvoxelbased broadbeam (nvbb) algorithm for intensity modulated radiation therapy dose calculation and plan optimization 
US7957507B2 (en)  20050228  20110607  Cadman Patrick F  Method and apparatus for modulating a radiation beam 
US20110202324A1 (en) *  20081020  20110818  Vanderbilt University  System and methods for accelerating simulation of radiation treatment 
US8229068B2 (en)  20050722  20120724  Tomotherapy Incorporated  System and method of detecting a breathing phase of a patient receiving radiation therapy 
WO2012099843A2 (en) *  20110117  20120726  Stc.Unm  System and methods for random parameter filtering 
US8232535B2 (en)  20050510  20120731  Tomotherapy Incorporated  System and method of treating a patient with radiation therapy 
US8363783B2 (en)  20070604  20130129  Oraya Therapeutics, Inc.  Method and device for ocular alignment and coupling of ocular structures 
US8442287B2 (en)  20050722  20130514  Tomotherapy Incorporated  Method and system for evaluating quality assurance criteria in delivery of a treatment plan 
US20130158879A1 (en) *  20060314  20130620  Xiaobo Sharon Hu  Methods and apparatus for hardware based radiation dose calculation 
US8506558B2 (en)  20080111  20130813  Oraya Therapeutics, Inc.  System and method for performing an ocular irradiation procedure 
US20140032185A1 (en) *  20100917  20140130  The University Of Texas M.D. Anderson Cancer Center  Gpubased fast dose calculator for cancer therapy 
US8663210B2 (en)  20090513  20140304  Novian Health, Inc.  Methods and apparatus for performing interstitial laser therapy and interstitial brachytherapy 
US8767917B2 (en)  20050722  20140701  Tomotherapy Incorpoated  System and method of delivering radiation therapy to a moving region of interest 
US20150199457A1 (en) *  20140114  20150716  Mitsubishi Electric Research Laboratories, Inc.  System and Method for Planning a Radiation Therapy Treatment 
US9443633B2 (en)  20130226  20160913  Accuray Incorporated  Electromagnetically actuated multileaf collimator 
US9731148B2 (en)  20050723  20170815  Tomotherapy Incorporated  Radiation therapy imaging and delivery utilizing coordinated motion of gantry and couch 
Families Citing this family (6)
Publication number  Priority date  Publication date  Assignee  Title 

WO2001085252A3 (en) *  20000505  20020523  Univ Washington  A method and apparatus for radiotherapy treatment planning 
US6714620B2 (en) *  20000922  20040330  Numerix, Llc  Radiation therapy treatment method 
US20080091388A1 (en) *  20030314  20080417  Failla Gregory A  Method for calculation radiation doses from acquired image data 
EP1901810B1 (en) *  20050526  20101027  Koninklijke Philips Electronics N.V.  Radiotherapeutic treatment planning incorporating functional imaging information 
US8125813B2 (en)  20050616  20120228  Best Medical International, Inc.  Variance reduction simulation system, program product, and related methods 
US7804936B2 (en) *  20070628  20100928  Siemens Medical Solutions Usa, Inc.  Doseguided radiation therapy using cone beam CT 
Citations (3)
Publication number  Priority date  Publication date  Assignee  Title 

WO1999040523A1 (en)  19980209  19990812  University Of Southampton  Treatment planning method and apparatus for radiation therapy 
US20020027971A1 (en) *  20000505  20020307  Deasy Joseph O.  Method and apparatus for radiotherapy treatment planning 
US20020106054A1 (en) *  20000922  20020808  Numerix, Llc.  Radiation therapy treatment method 
Patent Citations (4)
Publication number  Priority date  Publication date  Assignee  Title 

WO1999040523A1 (en)  19980209  19990812  University Of Southampton  Treatment planning method and apparatus for radiation therapy 
US6301329B1 (en)  19980209  20011009  The University Of Southampton  Treatment planning method and apparatus for radiation therapy 
US20020027971A1 (en) *  20000505  20020307  Deasy Joseph O.  Method and apparatus for radiotherapy treatment planning 
US20020106054A1 (en) *  20000922  20020808  Numerix, Llc.  Radiation therapy treatment method 
NonPatent Citations (6)
Title 

Buffa, F. et al. "Influence of Statistical Fluctuations in MonteCarlo Dose Calculations". Medical Physics, vol. 26, No. 6, Jun. 1999, pp. 11201120, XP001031982. 
Jiang, Steve B. et al. "Removing the effect of statistical uncertainty on dosevolume histograms from Monte Carlo dose calculations". Phys. Med. Biol., vol. 45, 2000, pp. 21512161. 
Mackie, T. Rockwell et al. "Tomotherapy: Optimized Planning and Delivery of Radiation Therapy". International Journal of Imaging Systems and Technology, Spring 1995, vol. 6, No. 1, pp. 4355, XP000620333. 
Petti, P.L. et al. "Investigation of buildup dose from electron contamination of clinical photon beams". Med. Phys., vol. 10, No. 1, Jan./Feb. 1983, pp. 1824. 
Sempau, J. et al. "Towards the elimination of Monte Carlo statistical fluctuation from dose volume histograms for radiotherapy treatment planning". Physics in Medicine and Biology, Jan. 2000, IOP Publishing, UK, vol. 45, No. 1, pp. 131157, XP002185777. 
Tabata, T; Andreo, P.; Ito, R. EnergyDeposition Distributions in Materials Irradiated by PlaneParallel Electron Beams with Energies Between 0.1 and 100 MeV. Atomic Data and Nuclear Data Tables, vol. 56, 1994 pp. 105131. 
Cited By (105)
Publication number  Priority date  Publication date  Assignee  Title 

US20050201516A1 (en) *  20020306  20050915  Ruchala Kenneth J.  Method for modification of radiotherapy treatment delivery 
US8406844B2 (en)  20020306  20130326  Tomotherapy Incorporated  Method for modification of radiotherapy treatment delivery 
US20040015073A1 (en) *  20020426  20040122  Michael Schell  Target repositioning error correction filter and method 
US20090063110A1 (en) *  20030314  20090305  Transpire,Inc.  Brachytherapy dose computation system and method 
US20060085223A1 (en) *  20041006  20060420  Jean Anderson  System and user interface for presenting treatment information 
US7970624B2 (en) *  20041006  20110628  Siemens Medical Solutions Usa, Inc.  System and user interface for presenting treatment information 
US8775207B2 (en) *  20050202  20140708  Siemens Medical Solutions Usa, Inc.  Integrated treatment planning and scheduling system 
US20060173725A1 (en) *  20050202  20060803  Michele Abraham  Integrated treatment planning and scheduling system 
US7957507B2 (en)  20050228  20110607  Cadman Patrick F  Method and apparatus for modulating a radiation beam 
US8232535B2 (en)  20050510  20120731  Tomotherapy Incorporated  System and method of treating a patient with radiation therapy 
WO2007014027A3 (en) *  20050722  20081002  Robert L Cravens  System and method of remotely analyzing operation of a radiation therapy system 
US8229068B2 (en)  20050722  20120724  Tomotherapy Incorporated  System and method of detecting a breathing phase of a patient receiving radiation therapy 
US20070104316A1 (en) *  20050722  20070510  Ruchala Kenneth J  System and method of recommending a location for radiation therapy treatment 
US20070041496A1 (en) *  20050722  20070222  Olivera Gustavo H  System and method of remotely analyzing operation of a radiation therapy system 
US20070041499A1 (en) *  20050722  20070222  Weiguo Lu  Method and system for evaluating quality assurance criteria in delivery of a treatment plan 
US20070041495A1 (en) *  20050722  20070222  Olivera Gustavo H  Method of and system for predicting dose delivery 
US7839972B2 (en)  20050722  20101123  Tomotherapy Incorporated  System and method of evaluating dose delivered by a radiation therapy system 
US8442287B2 (en)  20050722  20130514  Tomotherapy Incorporated  Method and system for evaluating quality assurance criteria in delivery of a treatment plan 
WO2007014027A2 (en) *  20050722  20070201  Tomotherapy Incorporated  System and method of remotely analyzing operation of a radiation therapy system 
US7773788B2 (en)  20050722  20100810  Tomotherapy Incorporated  Method and system for evaluating quality assurance criteria in delivery of a treatment plan 
US8767917B2 (en)  20050722  20140701  Tomotherapy Incorpoated  System and method of delivering radiation therapy to a moving region of interest 
US9731148B2 (en)  20050723  20170815  Tomotherapy Incorporated  Radiation therapy imaging and delivery utilizing coordinated motion of gantry and couch 
US20130158879A1 (en) *  20060314  20130620  Xiaobo Sharon Hu  Methods and apparatus for hardware based radiation dose calculation 
US8494115B2 (en) *  20060314  20130723  The University Of Notre Dame Du Lac  Methods and apparatus for hardware based radiation dose calculation 
US7680245B2 (en)  20061016  20100316  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US7693259B2 (en)  20061016  20100406  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US7693258B2 (en)  20061016  20100406  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US8761336B2 (en)  20061016  20140624  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US7697663B2 (en)  20061016  20100413  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US8837675B2 (en)  20061016  20140916  Oraya Therapeutics, Inc.  Ocular radiosurgery 
US7680244B2 (en)  20061016  20100316  Oraya Therapeutics, Inc.  Ocular radiosurgery 
US8855267B2 (en)  20061016  20141007  Oraya Therapeutics, Inc.  Orthovoltage radiosurgery 
US7564946B2 (en)  20061016  20090721  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US7822175B2 (en)  20061016  20101026  Oraya Therapeutics, Inc.  Portable orthovoltage radiotherapy 
US7535991B2 (en)  20061016  20090519  Oraya Therapeutics, Inc.  Portable orthovoltage radiotherapy 
US8442185B2 (en)  20061016  20130514  Oraya Therapeutics, Inc.  Orthovoltage radiosurgery 
US8995618B2 (en)  20061016  20150331  Oraya Therapeutics, Inc.  Portable orthovoltage radiotherapy 
US7912178B2 (en)  20061016  20110322  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US8320524B2 (en)  20061016  20121127  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US8189739B2 (en)  20061016  20120529  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US8073105B2 (en)  20061016  20111206  Oraya Therapeutics, Inc.  Ocular radiosurgery 
US20080187100A1 (en) *  20061016  20080807  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US20080089481A1 (en) *  20061016  20080417  Oraya Therapeutics, Inc.  Portable orthovoltage radiotherapy 
US20080089480A1 (en) *  20061016  20080417  Oraya Therapeutics, Inc.  Portable orthovoltage radiotherapy 
US7496174B2 (en)  20061016  20090224  Oraya Therapeutics, Inc.  Portable orthovoltage radiotherapy 
US8180021B2 (en)  20061016  20120515  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US8094779B2 (en)  20061016  20120110  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US8059784B2 (en)  20061016  20111115  Oraya Therapeutics, Inc.  Portable orthovoltage radiotherapy 
US8611497B2 (en)  20061016  20131217  Oraya Therapeutics, Inc.  Portable orthovoltage radiotherapy 
US8295437B2 (en)  20061213  20121023  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US7978819B2 (en)  20061213  20110712  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US8787524B2 (en)  20061213  20140722  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US7978818B2 (en)  20061213  20110712  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US8306186B2 (en)  20061213  20121106  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US7961845B2 (en)  20061213  20110614  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US8229069B2 (en)  20061213  20120724  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US8229073B2 (en)  20061213  20120724  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US9272161B2 (en)  20061213  20160301  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US7620147B2 (en)  20061213  20091117  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US8238517B2 (en)  20061213  20120807  Oraya Therapeutics, Inc.  Orthovoltage radiotherapy 
US7953203B2 (en)  20070409  20110531  Oraya Therapeutics, Inc.  Orthovoltage radiosurgery 
US8457277B2 (en)  20070409  20130604  Oraya Therapeutics, Inc.  Orthovoltage radiosurgery 
US7693260B2 (en)  20070409  20100406  Oraya Therapeutics, Inc.  Orthovoltage radiosurgery 
US7912179B2 (en)  20070409  20110322  Oraya Therapeutics, Inc.  Orthovoltage radiosurgery 
US8184772B2 (en)  20070409  20120522  Oraya Therapeutics, Inc.  Orthovoltage radiosurgery 
US8923479B2 (en)  20070604  20141230  Oraya Therapeutics, Inc.  Method and device for ocular alignment and coupling of ocular structures 
US8363783B2 (en)  20070604  20130129  Oraya Therapeutics, Inc.  Method and device for ocular alignment and coupling of ocular structures 
US8630388B2 (en)  20070604  20140114  Oraya Therapeutics, Inc.  Method and device for ocular alignment and coupling of ocular structures 
US20090110145A1 (en) *  20071025  20090430  Tomotherapy Incorporated  Method for adapting fractionation of a radiation therapy dose 
US8222616B2 (en)  20071025  20120717  Tomotherapy Incorporated  Method for adapting fractionation of a radiation therapy dose 
US8848869B2 (en)  20071223  20140930  Oraya Therapeutics, Inc.  Methods and devices for detecting, controlling, and predicting radiation delivery 
US8503609B2 (en)  20071223  20130806  Oraya Therapeutics, Inc.  Methods and devices for detecting, controlling, and predicting radiation delivery 
US7801271B2 (en)  20071223  20100921  Oraya Therapeutics, Inc.  Methods and devices for orthovoltage ocular radiotherapy and treatment planning 
US8494116B2 (en)  20071223  20130723  Oraya Therapeutics, Inc.  Methods and devices for orthovoltage ocular radiotherapy and treatment planning 
US7792249B2 (en)  20071223  20100907  Oraya Therapeutics, Inc.  Methods and devices for detecting, controlling, and predicting radiation delivery 
US9025727B2 (en)  20071223  20150505  Oraya Therapeutics, Inc.  Methods and devices for orthovoltage ocular radiotherapy and treatment planning 
US8920406B2 (en)  20080111  20141230  Oraya Therapeutics, Inc.  Device and assembly for positioning and stabilizing an eye 
US8506558B2 (en)  20080111  20130813  Oraya Therapeutics, Inc.  System and method for performing an ocular irradiation procedure 
US8512236B2 (en)  20080111  20130820  Oraya Therapeutics, Inc.  System and method for positioning and stabilizing an eye 
US20090234175A1 (en) *  20080311  20090917  Cynthia Maier  System and method for imageguided therapy planning and procedure 
US8333685B2 (en) *  20080311  20121218  Hologic, Inc.  System and method for imageguided therapy planning and procedure 
WO2009127747A1 (en)  20080414  20091022  Gmv Aerospace And Defence S.A.  Planning system for intraoperative radiation therapy and method for carrying out said planning 
US8880145B2 (en)  20080414  20141104  Gmv Aerospace And Defence S.A.  Planning system for intraoperative radiation therapy and method for carrying out said planning 
US20110052036A1 (en) *  20080414  20110303  Gmv Aerospace And Defence S.A.  Planning system for intraoperative radiation therapy and method for carrying out said planning 
US9592408B2 (en)  20080424  20170314  Koninklijke Philips N.V.  Dosevolume kernel generation 
US20110029329A1 (en) *  20080424  20110203  Koninklijke Philips Electronics N.V.  Dosevolume kernel generation 
US20100054413A1 (en) *  20080828  20100304  Tomotherapy Incorporated  System and method of calculating dose uncertainty 
US20100053208A1 (en) *  20080828  20100304  Tomotherapy Incorporated  System and method of contouring a target area 
US8803910B2 (en)  20080828  20140812  Tomotherapy Incorporated  System and method of contouring a target area 
US8913716B2 (en)  20080828  20141216  Tomotherapy Incorporated  System and method of calculating dose uncertainty 
US8363784B2 (en)  20080828  20130129  Tomotherapy Incorporated  System and method of calculating dose uncertainty 
US20110202324A1 (en) *  20081020  20110818  Vanderbilt University  System and methods for accelerating simulation of radiation treatment 
US9044601B2 (en) *  20081020  20150602  Dr. Fred J. Currell  System and methods for accelerating simulation of radiation treatment 
US20100228116A1 (en) *  20090303  20100909  Weiguo Lu  System and method of optimizing a heterogeneous radiation dose to be delivered to a patient 
US8663210B2 (en)  20090513  20140304  Novian Health, Inc.  Methods and apparatus for performing interstitial laser therapy and interstitial brachytherapy 
US20110122997A1 (en) *  20091030  20110526  Weiguo Lu  Nonvoxelbased broadbeam (nvbb) algorithm for intensity modulated radiation therapy dose calculation and plan optimization 
US8401148B2 (en)  20091030  20130319  Tomotherapy Incorporated  Nonvoxelbased broadbeam (NVBB) algorithm for intensity modulated radiation therapy dose calculation and plan optimization 
US9495513B2 (en) *  20100917  20161115  William Marsh Rice University  GPUbased fast dose calculator for cancer therapy 
US20140032185A1 (en) *  20100917  20140130  The University Of Texas M.D. Anderson Cancer Center  Gpubased fast dose calculator for cancer therapy 
US8712180B2 (en)  20110117  20140429  Stc.Unm  System and methods for random parameter filtering 
WO2012099843A2 (en) *  20110117  20120726  Stc.Unm  System and methods for random parameter filtering 
WO2012099843A3 (en) *  20110117  20121026  Stc.Unm  System and methods for random parameter filtering 
US9443633B2 (en)  20130226  20160913  Accuray Incorporated  Electromagnetically actuated multileaf collimator 
US9251302B2 (en) *  20140114  20160202  Mitsubishi Electric Research Laboratories, Inc.  System and method for planning a radiation therapy treatment 
US20150199457A1 (en) *  20140114  20150716  Mitsubishi Electric Research Laboratories, Inc.  System and Method for Planning a Radiation Therapy Treatment 
Also Published As
Publication number  Publication date  Type 

WO2001085252A2 (en)  20011115  application 
US20020027971A1 (en)  20020307  application 
WO2001085252A3 (en)  20020523  application 
Similar Documents
Publication  Publication Date  Title 

Wilson et al.  Noise properties of the EM algorithm. II. Monte Carlo simulations  
Nuyts et al.  Simultaneous maximum a posteriori reconstruction of attenuation and activity distributions from emission sinograms  
Hon et al.  A radial basis function method for solving options pricing models  
CRD  Statistical methods for tomographic image reconstruction  
Erdogan et al.  Monotonic algorithms for transmission tomography  
US6907102B1 (en)  Iterative reconstruction methods for multislice computed tomography  
Guan et al.  A projection access order for speedy convergence of ART (algebraic reconstruction technique): a multilevel scheme for computed tomography  
Vandenberghe et al.  Iterative reconstruction algorithms in nuclear medicine  
Silverman et al.  A smoothed EM approach to indirect estimation problems, with particular, reference to stereology and emission tomography  
Frey et al.  A practical method for incorporating scatter in a projectorbackprojector for accurate scatter compensation in SPECT  
Medoff et al.  Iterative convolution backprojection algorithms for image reconstruction from limited data  
Elbakri et al.  Statistical image reconstruction for polyenergetic Xray computed tomography  
Frey et al.  Application of taskbased measures of image quality to optimization and evaluation of threedimensional reconstructionbased compensation methods in myocardial perfusion SPECT  
Rantala et al.  Waveletbased reconstruction for limitedangle Xray tomography  
US6714620B2 (en)  Radiation therapy treatment method  
Li et al.  Reconstruction for proton computed tomography by tracing proton trajectories: A Monte Carlo study  
Kaufman  Maximum likelihood, least squares, and penalized least squares for PET  
US20060285640A1 (en)  Variance reduction simulation system, program product, and related methods  
US4991093A (en)  Method for producing tomographic images using direct Fourier inversion  
US5576548A (en)  Nuclear imaging enhancer  
Cho  General views on 3D image reconstruction and computerized transverse axial tomography  
Ye et al.  Nonlinear multigrid algorithms for Bayesian optical diffusion tomography  
Zanella et al.  Efficient gradient projection methods for edgepreserving removal of Poisson noise  
Wells et al.  Analytical calculation of photon distributions in SPECT projections  
US7203276B2 (en)  Xray scatter image reconstruction by balancing of discrepancies between detector responses, and apparatus therefor 
Legal Events
Date  Code  Title  Description 

AS  Assignment 
Owner name: WASHINGTON UNIVERSITY, MISSOURI Free format text: ASSIGNMENT OF ASSIGNORS INTEREST;ASSIGNORS:DEASY, JOSEPH O.;WICKERHAUSER, MLADEN VICTOR;REEL/FRAME:012127/0331;SIGNING DATES FROM 20010607 TO 20010620 

CC  Certificate of correction  
FPAY  Fee payment 
Year of fee payment: 4 

REMI  Maintenance fee reminder mailed  
LAPS  Lapse for failure to pay maintenance fees  
FP  Expired due to failure to pay maintenance fee 
Effective date: 20120914 