Method and apparatus for radiotherapy treatment planning

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

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

[0004] 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 non-tumorous tissue is minimized.

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

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

[0007] 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⁷), 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.

[0008] 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 well-defined 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.

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

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

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

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

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

[0014] 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:1765-1779), 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.

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

[0016] 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 non-parametric regression smoothing. The digital filtering techniques could employ Binomial/Gaussian filters or local-least squares filters.

[0017] 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 local-least squares filters. If employed, the Binomial/Gaussian filters could use a higher degree polynomial for less aggressive denoising. If employed, the local-least squares filters could use a smaller neighborhood for less aggressive denoising.

[0018] 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 batch-averaged wavelet denoising (BAWD) algorithm.

[0019] 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 pre-determined 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 pre-determined stopping criteria is reached.

[0020] 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 least-squares fit then extracting the fit to the center element of the submatrix as the least-squares fit, and repeating these acts for each element of the Monte Carlo data distribution D0 to complete the development of the denoised estimate.

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

[0022] In still another aspect of the present invention, developing the denoised estimate uses a statistically-based smoothing method. The statistically-based smoothing method could be the Loess method.

[0023] In yet another aspect of the present invention, developing the denoised estimate uses a curve fitting method employing Hermite polynomials for the least-squares 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 Fourier-based smoothing method. The Fourier-based smoothing method could employ a low-pass filter or a Wiener filter.

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

[0025] In a still further aspect of the present invention, the Monte Carlo data distribution D0 is regridded onto a non-cartesian grid and developing the denoised estimate occurs on the non-cartesian grid with the resulting denoised dose distribution D1 being transformed back to a cartesian grid. The non-cartesian grid could be a diverging fan beam grid.

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

[0027] In another aspect of the present invention, developing the denoised estimate uses a computational technique that reduces random noise and accurately estimates underlying signal.

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

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

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

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

[0032] FIG. 1 illustrates a digital step-heterogeneity 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;

[0033] FIG. 2 illustrates profile plots demonstrating digital filtering techniques to a uniform water phantom in accordance with the present invention;

[0034] FIG. 3 illustrates profile plots demonstrating digital filtering techniques to the step-heterogeneity phantom in accordance with the present invention;

[0035] FIG. 4 illustrates contour plots of dose distributions resulting from the uniform water phantom and the heterogeneity phantom in accordance with the present invention;

[0036] FIG. 5 illustrates application of spatially adaptive iterative filtering (SAIF) denoising algorithm to the heterogeneity phantom in accordance with the present invention;

[0037] FIG. 6 illustrates denoising of the heterogeneity phantom using wavelet denoising algorithm with 9,7 biorthogonal filters in accordance with the present invention;

[0038] FIG. 7 illustrates denoising of the heterogeneity phantom using wavelet shrinkage denoising with third-order ‘Coiflet’ wavelets and soft thresholding in accordance with the present invention; and

[0039] FIG. 8 illustrates wavelet denoising of three-dimensional photon dose distributions, in accordance with the present invention.

[0040] 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).

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

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

[0043] The Monte Carlo Method and Noise Characteristics

[0044] 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 square-root 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 square-root transform (i.e. taking the square-root of all dose values) is the appropriate variance stabilizing transform for such statistics (Poisson statistics). Hence the model:

[0045] Measured Dose Array^1/2=True Dose Array^1/2+Gaussian Noise Array  (1)

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

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

[0048] Denoising of Monte Carlo Electron Beam Dose Distributions via Digital Filtering

[0049] 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 two-dimensional filtering and image adaptive methods, although it will be apparent to those skilled in the art that the present invention is equally applicable to three-dimensional filtering and image adaptive methods. The filtering techniques are applied to a variance stabilized (square-root) dose image which is transformed back to dose-space (i.e. squared) after filtering.

[0050] Denoising is aimed at computing a new image h, having elements h_i,j, such that $\begin{array}{cc}{h}_{i,j}=\sum _{m=-k}^k\sum _{n=-k}^k{a}_{m,n}{g}_{i+n,j+m}& \left(2\right)\end{array}$

[0051] where the a_m,nvalues 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 (low-pass filters) and local-least squares (Savitzky-Golay) filters.

[0052] Binomial/Gaussian filters are constructed in 1-D 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.

[0053] In 2-D, 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 diffusion-like behavior and a Gaussian-like Fourier transfer functions. Hereafter, these filters are simply referred to as Gaussian filters.

[0054] Another interesting class of filters arises from the application of local least-squares 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: $\begin{array}{cc}h\left(x,y;a_0,\dots ,b_0,\dots ,b_n\right)=\sum _{p=0}^n{a}_p{x}^p\sum _{q=0}^m{b}_q{y}^q,& \left(3\right)\end{array}$

[0055] 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).

[0056] This type of filter is called a Local Least-Squares filter (LLS) and is also known as a Savitzky-Golay filter in applications to chemical spectral analysis. Previous applications of LLS filters to 2-D 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 3-D).

[0057] 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 step-heterogeneity phantom of FIG. 1, are through data column 60 (one of the two center-most columns in the 2-D 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 dip-and-rise is evident for all filters except LLS_W5_O3, where it cannot be visually detected. The blurring is worst for Gauss_W5_O4.

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

[0059] Denoising using digital filtering has the following advantages:

[0060] 1. Improved visualization of isodose lines (FIG. 3). The ability to visually inspect the dose distribution is significantly improved using denoising filtering.

[0061] 2. Suppression of spurious dose excursions (noise), as seen in the profile plots (FIGS. 2-3).

[0062] There is a trade-off 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.

[0063] Spatially Adaptive Iterative Filtering (SAIF)

[0064] In another aspect of the present invention, a one-parameter 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 chi-square comparison, over a local region, between the digital filter result and a ‘lightly’ denoised result which is known to be not too aggressive.

[0065] The method of spatially adaptive iterative filtering (SAIF):

[0066] Step 1. Begin with a safely denoised image (no significant bias), using the LLS_W5_O3 filter. The result is denoted F₀.

[0067] Step 2. Set the current best dose image (CBD) equal to F₀.

[0068] Step 3. Loop over increasingly aggressive filters:

[0069] Step 3.1. Denoise with the i'th filter, yielding F₁.

[0070] Step 3.2. Form a matrix of ‘local chi-square values,’ where each value is the chi-square comparison between the F₁doses and the F₀doses over a small region centered about that voxel.

[0071] Step 3.3. For each position in the chi-square image:

[0072] Step 3.3.1. If the local chi-square is below a preset threshold, replace that position in CBD with dose at that position in F₁.

[0073] FIG. 5 illustrates an example using SAIF. The local chi-square 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 chi-square 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 chi-square histogram shown in the second panel. The histogram shows how the local chi-square 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 well-approximated by the SAIF algorithm.

[0074] Wavelet Denoising

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

[0076] In one embodiment, a Wavelet Shrinkage Denoising Algorithm with 9,7-biorthogonal filters is used, the process summarized as follows:

[0077] Step 1. Apply the square-root variance stabilizing transformation.

[0078] Step 2. Spin-cycling: Create 8 extra matrices by shifting the original by one pixel value in all the nearest neighbor directions. (This eliminates grid orientation artifacts.)

[0079] Step 3. On the original plus the 8 shifted arrays do the following:

[0080] Step 3a. Transform: perform the forward wavelet transform using the 9,7-biorthogonal filters.

[0081] 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).

[0082] Step 3c. Inverse transform: perform the 9,7 inverse wavelet transform.

[0083] Step 4. Back-shift the 8 shifted arrays.

[0084] Step 5. Average the back-shifted arrays and the unshifted array.

[0085] Step 6. Inverse variance stabilization: square the resulting array values (step). This is the denoised dose distribution.

[0086] FIG. 6 illustrates denoising of the heterogeneity dose distribution of FIG. 1 using the wavelet denoising algorithm with 9,7-biorthogonal filters. Wavelet coefficient values less than 0.1% of the maximum wavelet coefficient value were set to zero. A computationally efficient lifting-scheme 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 low-density 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.

[0087] 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 step-heterogeneity phantom at a depth of 2.55 cm. Plots are through column 60 (one of the two center-most columns in the 2-D 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.

[0088] FIG. 8 illustrates wavelet denoising for 3-D photon dose distributions. The Pencil beam (PB) is denoised by denoising 2-D slices of the dose distribution in the x-z plane (i.e., a transverse-“beam direction” plane) using the 9,7-biorthogonal 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 (0-3 cm depth), aluminum (3-5 cm), homogeneous lung-like media (5-12 cm), and water (12-32 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 high-density central axis depth-dose with 1.2 billion source particles, the next shows a central axis depth-dose 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.

[0089] A thorough description of wavelet denoising using 9,7-biorthogonal 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).

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

[0091] 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,7-biorthogonal wavelet is chosen, as used by the draft JPEG 2000 image compression algorithm, since images are about as smooth as true underlying dose distributions.

[0092] Approximations of Dose Distributions

[0093] Monte Carlo radiotherapy simulations produce sampled dose distributions on a finite grid. A simulated two-dimensional dose distribution on an M×N grid may be represented by a non-negative 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 $\left(\frac{m}{M},\frac{n}{N}\right),$

[0094] given by the formula $\phi \left(x,y\right)=\sum _{m=0}^{M-1}\sum _{n=0}^{N-1}f\left(m,n\right){\phi }_{mn}\left(x,y\right).$

[0095] For example, characteristic functions of small regions could be taken: ${\phi }_{mn}^0\left(x,y\right)\stackrel{\mathrm{def}}{=}\{\begin{array}{cc}1,& \mathrm{if}\frac{m}{M}\le x<\frac{m+1}{M}\mathrm{and}\frac{n}{N}\le y<\frac{n+1}{N};\\ 0,& \mathrm{otherwise}.\end{array}$

[0096] A superposition of basic densities φ⁰_mnis constant on $\frac{1}{M}x\frac{1}{N}$

[0097] rectangles of the unit square. That is not a smooth or even continuous function.

[0098] Normal probability densities $N\left[\mu ,\sigma \right]\left(x\right)\stackrel{\mathrm{def}}{=}\frac{1}{\sqrt{2\pi \sigma }}\exp \left(-{\frac{1}{2}\left[\frac{x-\mu }{\sigma }\right]}^2\right),$

[0099] of suitable mean μ and variance σ, give smooth basic functions: ${\phi }_{mn}^{\infty }\left(x,y\right)\stackrel{\mathrm{def}}{=}N\left[\frac{m+\frac{1}{2}}{M},\frac{1}{M}\right]\left(x\right)\times N\left[\frac{n+\frac{1}{2}}{N},\frac{1}{N}\right]\left(y\right).$

[0100] A superposition of densities φ^∞_mnis 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.

[0101] The represented φ will be as smooth as its basic densities, regardless of f, so a selection between φ⁰_mnand φ^∞_mnis made. Basic densities with two continuous derivatives that are “almost” confined to the unit square are selected.

[0102] Efficient Coding of Fluctuations as Wavelets

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

[0104] 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,7-biorthogonal symmetric wavelets can be used for computing the scale and amplitude of fluctuations.

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

[0106] 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,7-biorthogonal symmetric wavelets, but is not necessary if other wavelet families are used as the basis of the denoising.

[0107] Discrete Wavelet Transforms

[0108] The wavelet transform is identical to that of JPEG 2000. For two-dimensional and three-dimensional arrays, the transform is applied separately along each dimension, so only the one-dimensional version is described. Let x=x(k), 0≦k<K be the array of K values to be transformed. Generate a new-array Fx=Fx(k), 0≦k<K, by the following rule: $\begin{array}{cc}Fx\left(k\right)=\{\begin{array}{cc}\sum _{n=-4}^{4}h\left(n\right)X\left(k-n\right),& \mathrm{if}k\mathrm{is}\mathrm{even};\\ \sum _{n=-4}^{2}g\left(n\right)X\left(k-1-n\right),& \mathrm{if}k\mathrm{is}\mathrm{odd}.\end{array}& \left(1\right)\end{array}$

[0109] Here h=h(n) and g=g(n) are the filter sequences defining the 9,7-biorthogonal 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 whole-sample symmetric reflection, as

[0110] defined by the formula below $X\left(k\right)=\{\begin{array}{cc}x\left(k\right),& \mathrm{if}0\le k<K;\\ x\left(2K-2-k\right),& \mathrm{if}K\le k<2K-1.\end{array}$

[0111] (Note: the above table and equation could use a higher resolution capture, but this is not absolutely necessary)

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

[0113] 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:

[0114] 1X(k)←X(k)+α[X(k−1)+X(k+1)], for all odd k in the range 0<k<K;

[0115] 2. X(k)←X(k)+β[X(k−1)+X(k+1)] for all even k in the range 0≦k<K;

[0116] 3. X(k)←X(k)+γ[X(k−1)+X(k+1)], for all odd k in the range 0<k<K;

[0117] 4. X(k)←X(k)+δ[X(k−1)+X(k+1)] for all even k in the range 0≦k<K;

[0118] 5. Fx(k)←X(k)/ζ, for all even k in the range 0≦k<K;

[0119] 6. Fx(k)←ζX(k), for all even k in the range 0≦k<K;

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

[0121] To invert F, use {−α,−β,−γ,−δ; 1/ζ} and apply the operations in reverse order:

[0122] 1. Y(k)←Y(k)/ζ, for all even k in the range 0≦k<K;

[0123] 2. Y(k)←Y(k)/ζ, for all odd k in the range 0<k<K;

[0124] 3. Y(k)←Y(k)−δ[Y(k−1)+Y(k+1)], for all even k in the range 0≦k<K;

[0125] 4. Y(k)←Y(k)−γ[Y(k−1)+Y(k+1)], for all odd k in the range 0<k<K;

[0126] 5. F¹y(k)←Y(k)−β[Y(k−1)+Y(k+1)], for all even k in the range 0≦k<K;

[0127] 6. F¹y(k)←Y(k)−α[Y(k−1)+Y(k+1)], for all odd k in the range 0<k<K.

[0128] The analog of Equation 1 of this section for F¹y is $\begin{array}{cc}{F}^{-1}y\left(n\right)=\sum _{k\in K_0\left(n\right)}{h}^{\prime }\left(2k-n\right)Y\left(2k\right)+\sum _{k\in K_1\left(n\right)}{g}^{\prime }\left(2k-n\right)Y\left(2k+1\right),& \left(2\right)\end{array}$

[0129] 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₀(n)={k:−4≦2k−n≦2} and K₁(n)={k:−4≦2k−n≦4}, for n=0,1, . . . , K−1.

[0130] Each even-indexed 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 odd-indexed elements Fx(1), Fx(3), . . . , etc.

[0131] To compute a one-dimensional wavelet transform on a signal x, apply F, retain the odd-indexed values, extract the even-indexed values to another array of about half the length, call that new array x, and repeat until x contains just a single number:

[0132] 1. Let j=1.

[0133] 2. Compute y(k)=Fx(k), k=0,1, . . . , K−1;

[0134] 3. Extract D_j(k)=y(2k+1) for all indices k with 0<2k+1<K.

[0135] 4. Replace x(k)=y(2k) for all indices k with 0≦2k<K.

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

[0137] 6. If K>1, then increment j←j+1 and go to step 2. Otherwise, stop.

[0138] At the termination of the above algorithm, the discrete wavelet transform of the original signal is stored in the arrays D₁, D₂, . . . , D_j, plus y(0). The values in D_jare the amplitudes of fluctuations of the signal at scales of 2^Jgrid 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.

[0139] 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₂, . . . , D₁), the concatenation of the fluctuation arrays ill the reverse order of their computation.

[0140] Approximation and De-Noising

[0141] The root mean square sum is used to measure the size of functions on grids M. $\Vert f\Vert \stackrel{\mathrm{def}}{=}\sqrt{\frac{1}{\left|M\right|}\sum _{m\in M}{\left|f\left(m\right)\right|}^2}.$

[0142] 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∥.

[0143] Define the thresholding function T=T_∈[x] for a function x on M by letting a=max{|x(m)|:m∈M} and then setting $T_{\epsilon }\left[x\right]\left(m\right)=\{\begin{array}{cc}0,& \mathrm{if}\left|x\left(m\right)\right|<\epsilon a,\\ x\left(m\right),& \mathrm{otherwise}.\end{array}$

[0144] Thresholding with small ε is a nonlinear transformation, but is close to the identity in the sense that $\frac{\left|x\left(m\right)-T_{\epsilon }\left[x\right]\left(m\right)\right|}{\Vert x\Vert }\le \epsilon ,$

[0145] for every m∈M.

[0146] 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(𝑓)₀being the single average coefficient, previously called y(0). The following algorithm is performed:

[0147] 1. Translate the function 𝑓 to 𝑓_o, so that 𝑓_o(0)=f(o).

[0148] 2. Compute the wavelet transform y_o←W(𝑓_o).

[0149] 3. Apply ε—thresholding z_o(m)←T(y_o(m)) for all m≠0.

[0150] 4. Compute the inverse wavelet transform x_o←W⁻¹(z_o).

[0151] 5. Translate the function x_oto x, so that x(o)=x_o(0).

[0152] 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 one-dimensional case, O={−1, 0, +1}, while in the two-dimensional case O is the nine-point 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.

[0153] 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 variance-stabilizing square-root operation and get all estimator for s from ${\left(E\left(\sqrt{d},\epsilon \right)\right)}^2.$

[0154] The threshold ∈ that minimizes the RMS error $\Vert {\left(E\left(\sqrt{d},\epsilon \right)\right)}^2-s\Vert$

[0155] 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 similar-size simulations.

[0156] Threshold Dependence on Run Length

[0157] 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 $\Vert {\left(E\left(\sqrt{d},\epsilon \right)\right)}^2-s\Vert$

[0158] taken over the R/r runs d of length r that comprise the run that produced s.

[0159] Spatially Adaptive Wavelet Denoising Algorithm (SAWD):

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

[0161] The Method of spatially adaptive wavelet denoising (SAWD):

[0162] 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₀.

[0163] Step 2. Set the current best dose image (CBD) equal to F₀.

[0164] Step. 3 Loop over increasing thresholds (for example: [0.1, 0.15, 0.2, 0.25]):

[0165] Step 3.1. Denoise with the wavelet denoising wden9_—7 software and the i'th threshold, yielding F₁.

[0166] Step 3.2. Form a matrix of ‘local chi-square values,’ where each value is the chi-square comparison between the F₁doses and the F₀doses over a small region centered about that voxel.

[0167] Step 3.3. For each position in the chi-square image:

[0168] Step 3.3.1. If the local chi-square is below a preset threshold, replace that position in CBD with the value in F_i.

[0169] Batch-averaged Wavelet Denoising:

[0170] In another aspect of the present invention, a Batch-averaged wavelet algorithm is employed to denoise the raw Monte Carlo dose distribution.

[0171] The Method of Batch-averaged wavelet denoising (BAWD):

[0172] 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):

[0173] Step 1. Apply the square-root variance stabilizing transformation.

[0174] Step 2. Spin-cycling: Create 8 extra matrices by shifting the original by one pixel value in all the nearest neighbor directions. (This eliminates grid orientation artifacts.)

[0175] Step 3. On the original plus the 8 shifted arrays do the following:

[0176] Step 3a. Transform: perform the forward wavelet transform using the 9,7-biorthogonal filters.

[0177] 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).

[0178] Step 3c. Inverse transform: perform the 9,7 inverse wavelet transform.

[0179] Step 4. Back-shift the 8 shifted arrays.

[0180] Step 5. Average the back-shifted arrays and the unshifted array.

[0181] Step 6. Inverse variance stabilization: square the resulting array values (step). This is the denoised dose distribution.

[0182] Step 7. End loop over batches (go to Step 0 if more batches are required).

[0183] Step 8. Return the average of the denoised batch solutions as the final answer.

[0184] 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 above-described 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.