US9251302B2 - System and method for planning a radiation therapy treatment - Google Patents
System and method for planning a radiation therapy treatment Download PDFInfo
- Publication number
- US9251302B2 US9251302B2 US14/154,725 US201414154725A US9251302B2 US 9251302 B2 US9251302 B2 US 9251302B2 US 201414154725 A US201414154725 A US 201414154725A US 9251302 B2 US9251302 B2 US 9251302B2
- Authority
- US
- United States
- Prior art keywords
- volume
- matrix
- radiation
- correcting
- beams
- 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.)
- Active, expires
Links
- 238000000034 method Methods 0.000 title claims abstract description 51
- 238000001959 radiotherapy Methods 0.000 title claims abstract description 22
- 230000005855 radiation Effects 0.000 claims abstract description 119
- 239000011159 matrix material Substances 0.000 claims abstract description 66
- 230000001186 cumulative effect Effects 0.000 claims abstract description 45
- 238000009826 distribution Methods 0.000 claims abstract description 38
- 230000001678 irradiating effect Effects 0.000 claims abstract description 9
- 239000013598 vector Substances 0.000 claims description 17
- 238000003892 spreading Methods 0.000 claims description 7
- 230000007480 spreading Effects 0.000 claims description 5
- 210000001519 tissue Anatomy 0.000 description 17
- 238000005457 optimization Methods 0.000 description 13
- 238000004364 calculation method Methods 0.000 description 12
- 238000000151 deposition Methods 0.000 description 12
- 230000008021 deposition Effects 0.000 description 11
- 238000012937 correction Methods 0.000 description 6
- 238000010586 diagram Methods 0.000 description 6
- 210000000988 bone and bone Anatomy 0.000 description 5
- 230000010339 dilation Effects 0.000 description 5
- 230000006870 function Effects 0.000 description 5
- 230000008859 change Effects 0.000 description 4
- 239000002245 particle Substances 0.000 description 4
- 230000002441 reversible effect Effects 0.000 description 4
- 238000003860 storage Methods 0.000 description 4
- 238000002591 computed tomography Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 3
- 238000000342 Monte Carlo simulation Methods 0.000 description 2
- 238000005284 basis set Methods 0.000 description 2
- 238000004891 communication Methods 0.000 description 2
- 238000007796 conventional method Methods 0.000 description 2
- 230000015654 memory Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000002829 reductive effect Effects 0.000 description 2
- 238000003325 tomography Methods 0.000 description 2
- 206010028980 Neoplasm Diseases 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 230000006978 adaptation Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 201000011510 cancer Diseases 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 238000005315 distribution function Methods 0.000 description 1
- 238000010894 electron beam technology Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000002710 external beam radiation therapy Methods 0.000 description 1
- 230000004907 flux Effects 0.000 description 1
- 230000005251 gamma ray Effects 0.000 description 1
- 230000001788 irregular Effects 0.000 description 1
- 238000012886 linear function Methods 0.000 description 1
- 239000004973 liquid crystal related substance Substances 0.000 description 1
- 238000002595 magnetic resonance imaging Methods 0.000 description 1
- 230000003211 malignant effect Effects 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 229940050561 matrix product Drugs 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 230000036961 partial effect Effects 0.000 description 1
- 238000002600 positron emission tomography Methods 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 230000002250 progressing effect Effects 0.000 description 1
- 238000007493 shaping process Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
- 238000012800 visualization Methods 0.000 description 1
Images
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
- A61N5/00—Radiation therapy
- A61N5/10—X-ray therapy; Gamma-ray therapy; Particle-irradiation therapy
- A61N5/103—Treatment planning systems
- A61N5/1031—Treatment planning systems using a specific method of dose optimization
-
- G06F17/5009—
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
- G16H20/00—ICT specially adapted for therapies or health-improving plans, e.g. for handling prescriptions, for steering therapy or for monitoring patient compliance
- G16H20/40—ICT specially adapted for therapies or health-improving plans, e.g. for handling prescriptions, for steering therapy or for monitoring patient compliance relating to mechanical, radiation or invasive therapies, e.g. surgery, laser therapy, dialysis or acupuncture
Definitions
- This invention relates generally to radiation therapy treatment planning, and particularly to determining a radiation dose for the radiation therapy.
- Radiation therapy is used to treat malignant tissue, such as cancer cells.
- the radiation can have an electromagnetic form, such as high-energy photons, or a particulate form, such as an electron, proton, neutron, or alpha particles.
- Dose determination generally includes two parts: a source model and a transport model.
- the source model provides the incident fluence, which is the flux of the radiation integrated over time.
- the transport model determines the dose that results from the incident fluence.
- the transport model is a performance bottleneck.
- the three main transport methods in an order of increasing accuracy or decreasing performance are pencil beam, superposition and convolution, and Monte Carlo simulation. Superposition and convolution is the most common method for determining the radiation dose for external beam radiation therapy.
- a tissue volume is irradiated by a large number of directed pencil beams of radiation at various depths.
- the treatment planning volume which represents a region within the patient, is partitioned into a rectangular 3D grid of voxels.
- To validate a treatment plan it is necessary to determine the irradiation pattern of tens of thousands of beams for millions of voxels in the volume.
- To optimize a dose plan it is necessary to determine the adjoint of this operation.
- the radiation pattern can be modeled as a spatial distribution of radiation deposited by each beam.
- the relative intensity of radiation delivered to the i th voxel by the j th beam is represented by a positive number A ij and the table of numbers for all such beams and voxels is called a fluence matrix A ⁇ R V ⁇ P , wherein V is a number of voxels and P is a number of pencil beams.
- V is a number of voxels and P is a number of pencil beams.
- Some dose calculation methods determine a cumulative voxel radiation dose applied during the treatment without explicitly constructing the fluence matrix.
- Some conventional methods determine the cumulative voxel radiation dose via Monte Carlo simulation or via clinically validated approximations. For example, one method calculate the contributions of the pencil-beams one at a time, e.g., by ray-tracing 3D dose distribution each pencil beam, and accumulate the result into a 3D array representing the cumulative doses to all voxels, see U.S. Pat. No. 8,325,878.
- Some state-of-the-art dose optimization techniques use a sparse approximation of the fluence matrix, in which a subset of fluence values above a predetermined threshold are stored, and the remaining small fluence values are neglected.
- a sparse enough matrix that can be stored in a computer memory disregards about 1-4% of the total radiation.
- the conventional methods are not fast enough for a real time dose calculation. Further, those methods are not suitable to determine the action of the adjoint A T of the fluence matrix, which can be necessary to determine the gradient of the error in various dose optimization techniques.
- Some embodiments of invention are based on a realization that there is no need to directly approximate or accumulate an accurate cumulative voxel radiation dose suitable for clinical use. Instead, it is advantageous to efficiently determine an initial approximation of the cumulative voxel radiation dose that is inaccurate and may not be suitable for t clinical use, and then subsequently correct the initial approximation with a set of independent operations having an adjoint.
- the approximation can be decomposed into independent operations to allow efficient calculation using multiple parallel processors. Because the correcting operations are also decomposable, the initial approximation can be subsequently corrected. Because each operation has a transpose, the transpose can be applied to the cumulative voxel radiation dose in a reverse order to obtain quantities needed for the dose optimization, such as the gradient of the dose error with respect to beam weights.
- various embodiments determine can increase the determination of the cumulative voxel radiation dose in one or two orders of magnitude.
- one embodiment of the invention discloses a method for planning a radiation therapy treatment of a body represented by a volume of voxels.
- the method includes determining a radiation dose matrix representing a spatial distribution of a radiation dose of beams of radiation irradiating the volume with homogeneous stopping power, wherein the beams are collimated and are shifted copies of each other at each depth of the volume, wherein the spatial distribution is determined by calculating cumulative beam-axial doses in a single matrix-matrix multiplication and redistributing the cumulative beam-axial doses to all voxels in the volume using a convolution; and applying a set of correcting operations to the radiation dose matrix to produce a cumulative voxel radiation dose of the volume, wherein each correcting operation is linear, independent from another correcting operation and has a transform.
- the steps of the method are performed by at least one processor.
- Another embodiment discloses a method for planning a radiation therapy treatment of a body represented by a volume of voxels including determining a radiation dose matrix representing a spatial distribution of a radiation dose of a regular 3D grid of beams irradiating the volume with homogeneous stopping power; and applying a set of correcting operations to the radiation dose matrix to produce a cumulative voxel radiation dose of the volume.
- the set of correcting operation includes one or combination of an operation for correcting an assumption of uniformity of lateral distribution of energy of each collimated beam, an operation for correcting an assumption of the homogeneous stopping power of the volume, an operation for correcting an assumption of a parallel distribution of the collimated beams, an operation for correcting an assumption of a multi-port propagation of the collimated beams, and an operation for correcting an assumption for a cumulative stopping power of media between a beam source and the patient, wherein each correcting operation is linear, independent from another correcting operation and has a transform.
- the steps of the method are performed by a parallel processor.
- Yet another embodiment discloses a radiation therapy system comprising a parallel processor for determining an approximation of a radiation dose matrix representing a spatial distribution of a radiation dose of beams of radiation irradiating a body represented by a volume of voxels, wherein the spatial distribution is determined by calculating cumulative beam-axial doses in a single matrix-matrix multiplication and redistributing the cumulative beam-axial doses to all voxels in the volume using a convolution, and for correcting the approximation using a set of linear correcting operations.
- FIG. 1 is a schematic diagram of a radiation therapy system 100 according to one embodiment of the invention.
- FIG. 2 is a schematic of a principle employed by some embodiments of the invention.
- FIG. 3 is a block diagram of a method for radiation therapy treatment planning using beams of radiation irradiating a volume of a patient according to some embodiments of the invention
- FIG. 4 is a block diagram of a method for determining an approximation of the spatial distribution of the radiation dose according to one embodiment
- FIG. 5 is a schematic showing a set of correcting operations used by some embodiment of the invention to correct inaccuracy of initial approximation
- FIG. 6 is a schematic of a realization used by one embodiment of the invention.
- FIG. 7 is a schematic of a method for correcting the assumption of a parallel distribution of the collimated beams according to one embodiment of the invention.
- FIGS. 8A and 8B are flow charts for an exemplar method for determining a cumulative voxel radiation dose for radiation therapy treatment of the volume of the patient according to one embodiment of the invention.
- FIG. 1 is a schematic diagram of a radiation therapy system 100 according to one embodiment of the invention.
- the radiation therapy system 100 includes a radiation treatment planning system 101 , which further includes a parallel processor 102 .
- the parallel processor is adapted to receive input information concerning a body 105 having an intended radiation treatment volume that can be represented as a volume of voxels.
- the parallel processor 102 is also adapted to generate output information for providing radiation treatment to the intended radiation treatment volume of the body.
- the radiation treatment planning system 101 can further include a storage 107 , a display 108 , and input/output (I/O) devices and interfaces 109 .
- the storage 107 may be, for example, a hard disk drive, a CD-ROM drive, a DVD drive, a flash drive, etc.
- the display 108 may be, for example, a liquid crystal display (LCD), a cathode ray tube (CRT) monitor, a plasma display, etc.
- I/O device 109 may include, for example, a mouse, a keyboard, an interface for data transfer over a network or a data bus.
- the radiation therapy system 100 can further include a radiation treatment system 103 that is in communication with the radiation treatment planning system 101 .
- the radiation treatment system 103 can include a radiation source 106 that emits a directed beam of radiation for treatment.
- radiation sources may include, a X-ray source, a gamma ray source, an electron beam source, etc.
- the radiation source 106 may further comprise a multi-leaf collimator (MLC) to shape the beam. By adjusting the position of the leaves of the MLC, a radiotherapist can match the radiation field to a shape of the treatment volume of body. Other beam shaping and/or contouring can be included in some embodiments.
- the radiation source 106 can have a corresponding source model.
- the radiation system 103 may be controlled by the radiation treatment planning system 101 , for example, to deliver intensity modulated radiation energy and to conform radiation treatment to the shape of the intended radiation treatment volume.
- the radiation therapy system 100 can also include a diagnostic system 104 , in communication with the radiation treatment planning system 101 that generates empirical data of body 105 , e.g., a body of a human.
- the empirical data may be used as input information to the radiation treatment planning system 101 and the parallel processor 102 and may be used to determine the dose of radiation.
- the diagnostic system 104 can include sensors to obtain the empirical data of body 105 .
- An example diagnostic system may be a computed tomography (CT) scanner, a magnetic resonance imaging (MRI) scanner, a positron emission tomography (PET) scanner.
- CT computed tomography
- MRI magnetic resonance imaging
- PET positron emission tomography
- Some embodiments of the invention are based on a recognition that for particle beam dose calculation and optimization, it is necessary to rapidly determine cumulative voxel radiation dose using a matrix-vector products Aw and/or A T ⁇ , without explicitly constructing the matrix A, wherein T is the transpose operator, w is a vector of nonnegative beam weights and ⁇ is a vector of errors in the doses of the voxels.
- FIG. 2 shows a schematic of a principle employed by some embodiments of the invention.
- a volume 230 of a body of a patient is irradiated 220 by a large number of pencil beams, e.g., the beams 240 - 241 , that are directed at different points in the patient, and penetrate to different depths.
- the beams can be considered as translated copies of each other, and the 3D irradiation field of each beam can be modeled 240 as a trilinear outer product 250 of three vectors.
- the trilinear outer product can be used to rapidly approximate a linear superposition of many parallel beams directed at a rectangular grid of targets.
- some embodiments determine the trilinear outer product by calculating cumulative beam-axial doses in a single matrix-matrix multiplication and redistributing the doses to all voxels in the volume of the patient in a convolution, as described in more details below. Such calculations further increase the efficiency of the computation.
- the beams of radiation aimed at different depths have different radiation profiles.
- the irradiation pattern of each individual beam depends on the stopping and scattering power of the tissues the beam passes through.
- the beams usually radiate from a common source 210 and propagation of each beam increases nonlinearly.
- Some embodiments of the invention apply a sequence of corrections 260 to the dose approximation to account for all of these phenomena.
- Each of these correction operations has a transpose.
- some embodiments obtain quantities needed for dose optimization, such as the gradient of the dose error with respect to the beam weights.
- FIG. 3 shows a block diagram of a method for radiation therapy treatment planning using beams of radiation irradiating a volume of a patient.
- the method determines a cumulative voxel radiation dose for the volume, and each step of the method can be implemented by the parallel processor 102 .
- the method determines 310 a radiation dose matrix representing a spatial distribution of a radiation dose of a regular 3D grid of beams of radiation.
- the radiation dose matrix is approximated 315 based on the assumption that the beams are collimated, and collimated beams at each depth are shifted copies of each other for radiating the volume of the patient with homogeneous stopping power. Such an assumption allows for rapid determination of the approximation of the radiation dose on the parallel processor.
- the method applies 320 a set of correcting operations to the radiation dose matrix to produce a cumulative voxel radiation dose 325 of the volume.
- Each correcting operation is linear, independent from another correcting operation, and has a transpose.
- Some embodiments approximate a spatial distribution of a radiation dose based on assumption that all beams are of a given energy and are grid-translated copies of a single beam. For each beam, the density of energy deposition by the beam (as a function of depth z) is presumed to follow some energy profile, e.g., a Bragg profile B(z).
- FIG. 4 shows a block diagram of a method for determining the approximation of the spatial distribution of the radiation dose according to one embodiment.
- the spatial distribution is determined by calculating cumulative beam-axial doses in a single matrix-matrix multiplication 410 and redistributing the doses to all voxels in the volume of the patient in a convolution 420 .
- N be a column vector containing the values of the lateral distribution (e.g., a 1D Gaussian density function) and let N be a matrix in which each column is a spatially shifted copy of N, the shift matching the lateral spacing between beams.
- B be a vector containing the values of a Bragg profile and let B be a matrix of spatially shifted Bragg profiles (or let each column in B contain a unique energy profile for all beams in one energy layer).
- vec (UXV) (V T U) vec (X)
- Aw vec( B T mat( w )( N N )) where mat (•) reshapes its argument into a matrix that conforms to its left and right multiplicands.
- This enables is a very efficient mechanism to obtain the dose deposition Aw by multiplying 430 a vector of beam weights by a matrix of energy profiles, e.g., as a product B T mat (w).
- Each value in B T mat (w) is the cumulative radiation deposited into a specific Z plane by all beams that center on a specific X,Y coordinate.
- one embodiment substitute 440 values in voxels in the volume with corresponding values of the total radiation, i.e., arrange these values at their XYZ coordinates in a 3D array otherwise filled with zeroes, and then convolves 420 the volume with a 2D (XY) lateral spreading kernel.
- a 2D convolution can be replaced with two 1D convolutions.
- the convolving includes a first convolution with a first vector in a first lateral direction, followed by a second convolution with a second vector in a second lateral direction orthogonal to the first lateral direction, e.g., one convolution is in X direction followed by the convolution in Y direction.
- Both the matrix-vector product and the convolutions can be natively parallelized by a graphics processor unit (GPU).
- GPU graphics processor unit
- the approximation can be decomposed into independent operations to allow efficient calculation using multiple processors. Because the correcting operations are also decomposable, the initial approximation can be subsequently corrected in parallel with the help of multiple processors. Because each operation has a transpose, the transpose can be applied to the cumulative voxel radiation dose in reverse order to obtain quantities needed for the dose optimization, such as the gradient of the dose error with respect to the beam weights. However, the trilinear approximation of the cumulative voxel radiation dose is inaccurate based on one or combination of the following.
- the approximation assumes that transverse slices of all beams have identical, e.g., Gaussian, lateral distributions at all depths.
- the variance of the lateral distribution increases nonlinearly with increasing depth in the patient, so that beam has a unique dilation at each depth.
- the approximation assumes irradiation of the volume with a homogeneous stopping power, e.g., water.
- a homogeneous stopping power e.g., water.
- the tissues in the volume of the patient are inhomogeneous, i.e., voxels representing bones have a relatively high radiation stopping power, while voxels representing air have very little stopping power.
- Some embodiments of the invention are based on a realization that the inaccuracy of the approximation due to the variance of the stopping power in the volume can be reduced by stretching and contracting the voxels in the original homogeneous volume along the line of the propagation of the beam until the length of each voxel is inversely proportioned to its stopping power. For example, if the approximation of the cumulative voxel radiation dose is determined for the volume of water, a bone in the volume effectively shortens the beam, while the air effectively lengthens the beam. With such stretching and contracting operations, the dose deposition at the usual linear rate delivers the correct dose to each voxel.
- embodiments based on this realization account for inhomogeneous tissues by making corrections to the cumulative dose deposited in a homogeneous volume based on a redistribution matrix R that linearly maps irradiation values from those deposited into regular grid into those that would be deposited into a stretched grid.
- the redistribution matrix can be predetermined, e.g., based on Hounsfield radiodensity values associated with computer-aided tomography density readings of the patient.
- the application of the redistribution matrix can also be performed in parallel using multiple processors and has an adjoint operation.
- the approximation assumes irradiation of the volume with parallel beams.
- the beams can radiate from a common focal point.
- Some embodiments are based on a recognition that a Euclidean space can be warped so that all lines parallel to one axis meet at a focal point.
- some embodiments counter-warp the patient geometry and tissue density data.
- dose deposition and redistribution proceeds as described above.
- the planning volume is warped back to its original geometry.
- the warp and counter-warp are linear operators that can be handled by parallel computing operations or represented by a very sparse matrix that can be combined (or sequenced) with the rotation and heterogeneous tissue redistribution matrices.
- the warp and counter-warp operations change the volume of any subvolume but do not change the stopping power or cumulative radiation in that subvolume.
- the counter-warp can also be applied to the beam-spreading profiles to keep the geometry of the beams consistent with the geometry of the treatment volume.
- FIG. 5 shows a schematic of a set of correcting operations 501 used by some embodiment of the invention to correct inaccuracy of the approximation 315 .
- the set of correcting operations 510 includes one or combination of an operation for correcting an assumption of uniformity of lateral distribution of energy of each collimated beam 505 , an operation for correcting an assumption of the homogeneous stopping power of the volume 510 , an operation for correcting an assumption of a parallel distribution of the collimated beams 520 , an operation for correcting an assumption of a multi-port propagation of the collimated beams 530 , and an operation for correcting an assumption for a cumulative stopping power of media between a beam nozzle and the patient volume 540 .
- One embodiment approximates the exact sum (eq. 3.1) with a weighted superposition of a small number of differently varianced and re-weighted outer products.
- This embodiment is based on a realization that for each basis, there is a function giving modified weights u ij (w i B i (z), ⁇ i 2 (z), ⁇ j 2 (z)) and, thus, the exact sum can be approximated as ⁇ j J ( ⁇ i u ij ( w i B i ( z ), ⁇ i 2 ( z ), ⁇ j 2 ( z ))) ⁇ N (( x,y ), ⁇ j 2 ( z )). (3.3)
- the new weighting function u ij (•) can be obtained by solving a small and simple linear system of equations, as described below.
- the embodiment determines the approximation to be most accurate near the beam core, and weights the squared error by another 2D Gaussian field N((0,0), ⁇ 2 ).
- the approximation weights depend only on the depth and the superposition further simplifies to a form that is close to the original trilinear approximation in eqn. 3.2, indicating that dilation across all beams can be modeled by summing a small number of trilinear approximations, each with different transverse spreads: ⁇ j J ( ⁇ i w i v j ( z ) B i ( z )) ⁇ N (( x,y ), ⁇ j 2 ). (3.6)
- a weighted sum of 20 Gaussian fields with variances drawn randomly from the interval ⁇ i (z) ⁇ [3,6] mm can be approximated using a basis set of 4 Gaussian fields with fixed variances such that the maximum absolute error is less than 1/1000th of the exact peak value.
- the approximation can be made tighter by using more fields in the basis, or faster by using less.
- this correction operation uses a linear operator, i.e., its transpose is also available for dose optimization.
- voxels representing bone have a relatively high radiation stopping power, while voxels representing air have very little stopping power. To a first order, this can be modeled by stretching and contracting the voxels along the beam until they are inversely proportioned to their stopping power. Then, the dose deposition at the usual linear rate determines the correct dose for each voxel. Bone effectively shortens the beam; air effectively lengthens the beam.
- This “dose-volume equivalent” approximation usually requires serial ray-tracing of each beam progressing through the volume. The ray tracing can be accelerated using GPUs and reverse ray tracing techniques, but another two orders of magnitude speed-up is needed for real-time dose optimization. Furthermore, ray-tracing methods do not lend themselves to the gradient computation needed for dose optimization.
- FIG. 6 shows an illustration of a realization used by one embodiment of the invention.
- water-equivalent dosing 610 denser voxels are stretched to yield an irregular grid on which then radiation dosage is calculated.
- redistribution dosing 612 employed by the embodiment the radiation dosage is calculated on a regular grid and then redistributed to match the dosage pattern of the stretched voxels.
- one embodiment uses a perturbational procedure that is amenable to massively parallel computation.
- the embodiment precalculates a redistribution matrix R that linearly maps irradiation values from those deposited into regular grid into those that would be deposited into a stretched grid.
- This redistribution matrix corrects a dose deposited in a column of water to match the dose that would be deposited by ray-tracing through a column of inhomogeneous tissue.
- the (nonzero) values of R can be derived, e.g., from the Hounsfield radiodensity values associated with computer-aided tomography density readings.
- FIG. 0 illustrates.
- the dose deposition A and redistribution R are linear and can be combined to give the correct heterogeneous-tissue dose RAw.
- the redistribution is many-to-many, i.e., the irradiation values from several contiguous voxels in the initial dose calculation Aw grid may be linearly combined to obtain the total radiation absorbed by single voxel of bone or very dense tissue the heterogeneity-corrected dose RAw.
- irradiation in a single voxel in Aw can be distributed to many low-density voxels in RAw.
- the redistribution matrix R is sparse and highly compressible, i.e., for each voxel v i in RAw, one embodiment determines the shallowest and deepest voxels contributing voxels in Aw and the fractions of their irradiations contributed to v i . Consequently, in one embodiment, the redistribution matrix R can be represented with only four nonzero values per row providing an efficient storage and use by the GPU. For other kinds of radiation sources, the redistribution matrix may not be as sparse. In those cases, some embodiments express the redistribution as a piecewise linear function, such that the texture-mapping functions of the GPU can be utilized.
- the mathematically idealized dose calculation assumes that all beams are parallel; in reality, all beams radiate from a common focal point. To handle this, some embodiments based on a recognition that a Euclidean space can always be warped so that all lines parallel to one axis meet at the focal point. Similarly, there is a counter-warp that inverts this warp (except at the focal point, which is a singularity).
- FIG. 7 shows a schematic of a method for correcting the assumption of a parallel distribution of the collimated beams according to one embodiment of the invention.
- the embodiment counter-warp 710 the patient geometry 720 and tissue density data. Then dose deposition and redistribution proceeds as above.
- the counter-warping transforms 730 the diverging beam 740 into the parallel beams 750 used to determine the approximation 315 .
- the counter-warping is reversed after the radiation dose is determined and the planning volume is warped back to its original geometry.
- the warp and counter-warp operation change the volume of any subvolume of the volume but do not change the stopping power or cumulative radiation in that volume.
- the counter-warp can also be applied to the beam-spreading profiles to keep the geometry of the beams consistent with the geometry of the volume.
- the warp and counter-warp are linear operators that can be handled by native GPU operations or represented by a very sparse matrix that can be combined (or sequenced) with the rotation and heterogeneous tissue redistribution matrices.
- the same technique can be used to calculate the dose from beams that are not aligned with the Z-axis of the treatment planning volume, as is the case for multiport planning.
- the beams from the first port travel parallel to the Z axis of the treatment planning volume, and the beams of the second port travel along some other axis specified by a set of rotations, for example Euler angles (azimuth, elevation, twist).
- some embodiments predetermine a rotation matrix that redistributes the dose in each voxel of the original axis-aligned volume to the appropriate voxels in a counter-rotated volume.
- one embodiment overlies the original volume and the counter-rotated volume, such that for each voxel in the original volume, its fraction is overlap with each voxel in the counter-rotated volume.
- This approach produces another sparse matrix P with a small number of non-zeros per row that is easily handled by the GPU. In fact, this is exactly how GPUs do rotation in volume visualization. Because rotation is linear, the rotation can be combined with the redistribution operator for heterogeneous tissue.
- an independent rotation P k , redistribution R k , beam basis A k , and set of beam weights w k are associated with each port, yielding a cumulative all-ports dose of ⁇ k P k R k A k w k .
- the Lagrange multipliers are shared between all ports, because all ports cooperate to satisfy the constraints.
- the planning volume can represent the interior of the patient, but the beams are attenuated by the air, skin, tissues, and other media that lie between the planning volume and the beam source.
- the effect of these media on the beams can be modeled in dose optimization through the dose redistribution matrix R.
- R dose redistribution matrix
- one embodiment adds a special layer of voxels between the TPV and the beam nozzle.
- the dose-volume equivalent stopping power of each voxel in this layer is set to the cumulative dose-volume stopping equivalent of all the media on the line between that voxel and the nozzle.
- the R matrix is formed as before. If there are no constraints on these special voxels, the corresponding rows of the R matrix can be dropped in the optimization.
- FIGS. 8A and 8B show a flow chart of an exemplar method for determining a cumulative voxel radiation dose for radiation therapy treatment of the volume of the patient.
- the method tests 805 if the data representing density of the volume of the patient needs to be rotated. For example, if the radiation treatment system 103 includes multiple source of power, the data is counter-rotated 810 . Similarly, if the pencil beams of radiation are diverging 815 , the density data are counter-wrapped 820 as described above.
- the radiation deposited at each depth by each set of coaxial beams is determined 825 by multiplying the beam profile matrix by a matrix of beam weights. These deposition values are seeded 830 into the voxels along each beam line in the 3D voxel array. The rest of the 3D voxel array is zeroed 835 and the axial deposition values are spread throughout the volume 840 via planar convolutions with the lateral spread kernels. If the beams are known to dilate 850 , steps 825 - 840 are repeated with slightly different beam profile matrices 855 and the results are linearly combined 860 to finely approximate the true radiation deposition pattern.
- the deposited values are redistributed up and down the 3D volume 870 to model the varied stopping power of heterogeneous tissue.
- any compensation for diverging beams 875 is inverted 880 and similarly any compensation for patient rotation 885 is inverted 890 to restore the patient volume to its original geometry.
- the embodiments can be implemented in any of numerous ways.
- the embodiments may be implemented using hardware, software or a combination thereof.
- the software code can be executed on any suitable processor or collection of processors, whether provided in a single computer or distributed among multiple computers.
- processors may be implemented as integrated circuits, with one or more processors in an integrated circuit component.
- a processor may be implemented using circuitry in any suitable format.
- the processor can be connected to memory, transceiver, and input/output interfaces as known in the art.
- the various methods or processes outlined herein may be coded as software that is executable on one or more processors that employ any one of a variety of operating systems or platforms.
- the invention may be embodied as a computer readable medium other than a computer-readable storage medium, such as signals.
- program or “software” are used herein in a generic sense to refer to any type of computer code or set of computer-executable instructions that can be employed to program a computer or other processor to implement various aspects of the present invention as discussed above.
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Biomedical Technology (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Physics (AREA)
- General Engineering & Computer Science (AREA)
- Radiology & Medical Imaging (AREA)
- Life Sciences & Earth Sciences (AREA)
- Software Systems (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Veterinary Medicine (AREA)
- Public Health (AREA)
- General Health & Medical Sciences (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Pathology (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Animal Behavior & Ethology (AREA)
- Geometry (AREA)
- Evolutionary Computation (AREA)
- Computer Hardware Design (AREA)
- Radiation-Therapy Devices (AREA)
Abstract
Description
E(x 0 ,y 0 ,z 0,Δ)=∫z
E sep(x 0 ,y 0 ,z 0,Δ)=∫z
Aw=vec(B Tmat(w)(N N))
where mat (•) reshapes its argument into a matrix that conforms to its left and right multiplicands.
This enables is a very efficient mechanism to obtain the dose deposition Aw by multiplying 430 a vector of beam weights by a matrix of energy profiles, e.g., as a product BT mat (w). Each value in BT mat (w) is the cumulative radiation deposited into a specific Z plane by all beams that center on a specific X,Y coordinate. To calculate the dose, one
Σi w i B i(z)N((x,y),σi 2(z)), (3.1)
where wi is the weight on the ith beam, Bi(z) is the integral of the Bragg profile over the interval spanned by volume layer z, and N((x,y), σi 2(z)) is a 2D Gaussian field of variance σi 2(z) centered at (x,y). Calculating and summing all these fields is very expensive. In the trilinear outer product procedure described above, all variances are presumed equal (∀i,zσi 2(z)=σ2), and the sum is factored as
(Σi w i B i(z))·N((x,y),σ2). (3.2)
Σj J(Σi u ij(w i B i(z),σi 2(z),σj 2(z)))·N((x,y),σj 2(z)). (3.3)
E=∫∫N((x,y);(0,0),α2)·(N((x,y);(0,0),β2)−Σj J v j N((x,y);(0,0),σj 2))2 dxdy.
where each sum includes k=j.
Σj J{(Σi v ij(σi 2(z),σj 2(z))·w i B i(z))·N((x,y),σj 2(z)),},
where each term in braces can be determined by the original outer product procedure. The approximation weights vij can be pre-calculated (for every beam and z layer) and re-used for different beam weights wi.
Σj J(Σi w i v j(z)B i(z))·N((x,y),σj 2). (3.6)
Σj J(Σi w i {hacek over (B)} i(z))·N((x,y),σj 2). (3.7)
Claims (12)
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
US14/154,725 US9251302B2 (en) | 2014-01-14 | 2014-01-14 | System and method for planning a radiation therapy treatment |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
US14/154,725 US9251302B2 (en) | 2014-01-14 | 2014-01-14 | System and method for planning a radiation therapy treatment |
Publications (2)
Publication Number | Publication Date |
---|---|
US20150199457A1 US20150199457A1 (en) | 2015-07-16 |
US9251302B2 true US9251302B2 (en) | 2016-02-02 |
Family
ID=53521596
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
US14/154,725 Active 2034-08-01 US9251302B2 (en) | 2014-01-14 | 2014-01-14 | System and method for planning a radiation therapy treatment |
Country Status (1)
Country | Link |
---|---|
US (1) | US9251302B2 (en) |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US11147985B2 (en) * | 2017-02-02 | 2021-10-19 | Koninklijke Philips N.V. | Warm start initialization for external beam radiotherapy plan optimization |
US11541252B2 (en) * | 2020-06-23 | 2023-01-03 | Varian Medical Systems, Inc. | Defining dose rate for pencil beam scanning |
Families Citing this family (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US9251302B2 (en) * | 2014-01-14 | 2016-02-02 | Mitsubishi Electric Research Laboratories, Inc. | System and method for planning a radiation therapy treatment |
US10500414B2 (en) * | 2014-10-24 | 2019-12-10 | National Institutes For Quantum And Radiological Science And Technology | Irradiation planning apparatus and irradiation plan correction method |
US10417390B2 (en) * | 2015-06-30 | 2019-09-17 | Varian Medical Systems, Inc. | Methods and systems for radiotherapy treatment planning |
CN111494815B (en) * | 2020-05-14 | 2022-04-29 | 安徽慧软科技有限公司 | Three-dimensional dose calculation method, device and medium based on mixed variable-scale model |
Citations (16)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20020027971A1 (en) * | 2000-05-05 | 2002-03-07 | Deasy Joseph O. | Method and apparatus for radiotherapy treatment planning |
US6714620B2 (en) * | 2000-09-22 | 2004-03-30 | Numerix, Llc | Radiation therapy treatment method |
US20060259282A1 (en) * | 2003-03-14 | 2006-11-16 | Failla Gregory A | Deterministic computation of radiation transport for radiotherapy dose calculations and scatter correction for image reconstruction |
US7496173B2 (en) * | 2003-12-12 | 2009-02-24 | University Of Western Ontario | Method and system for optimizing dose delivery of radiation |
US7529339B2 (en) * | 2003-12-12 | 2009-05-05 | University Of Western Ontario | Method and system for optimizing dose delivery of radiation |
US7804935B2 (en) * | 2003-08-28 | 2010-09-28 | Henry Ford Health System | Fuzzy logic guided inverse treatment planning |
US7945022B2 (en) * | 2008-03-12 | 2011-05-17 | Sun Nuclear Corp. | Radiation therapy plan dose perturbation system and method |
US20120041685A1 (en) | 2010-07-13 | 2012-02-16 | Vanderbilt University | System and method for estimating radiation dose and distribution using medium-dependent-correction based algorithms (mdc) |
US8325878B2 (en) * | 2008-05-08 | 2012-12-04 | The Johns Hopkins University | Real-time dose computation for radiation therapy using graphics processing unit acceleration of the convolution/superposition dose computation method |
US8403821B2 (en) * | 2008-03-17 | 2013-03-26 | Mitsubishi Heavy Industries, Ltd. | Radiotherapy apparatus controller and radiation irradiating method |
US8422631B2 (en) * | 2007-12-07 | 2013-04-16 | Mitsubishi Heavy Industries, Ltd. | Radiation therapy planning apparatus and radiation therapy planning method |
US20130102830A1 (en) * | 2010-06-22 | 2013-04-25 | Karl Otto | System and method for estimating and manipulating estimated radiation dose |
US8682414B2 (en) * | 2009-11-27 | 2014-03-25 | Mitsubishi Heavy Industries, Ltd. | Radiation therapy apparatus control method and radiation therapy apparatus controller |
US20150038762A1 (en) * | 2011-12-13 | 2015-02-05 | Gsi Helmholtzzentrum Fuer Schwerionenforschung Gmbh | Method for preparing an irradiation plan for a moving target volume without a movement compensation |
US20150199457A1 (en) * | 2014-01-14 | 2015-07-16 | Mitsubishi Electric Research Laboratories, Inc. | System and Method for Planning a Radiation Therapy Treatment |
US20150202464A1 (en) * | 2014-01-23 | 2015-07-23 | Mitsubis | Multi-Criteria Optimization in Particle Beam Dose Optimization |
-
2014
- 2014-01-14 US US14/154,725 patent/US9251302B2/en active Active
Patent Citations (17)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20020027971A1 (en) * | 2000-05-05 | 2002-03-07 | Deasy Joseph O. | Method and apparatus for radiotherapy treatment planning |
US6792073B2 (en) * | 2000-05-05 | 2004-09-14 | Washington University | Method and apparatus for radiotherapy treatment planning |
US6714620B2 (en) * | 2000-09-22 | 2004-03-30 | Numerix, Llc | Radiation therapy treatment method |
US20060259282A1 (en) * | 2003-03-14 | 2006-11-16 | Failla Gregory A | Deterministic computation of radiation transport for radiotherapy dose calculations and scatter correction for image reconstruction |
US7804935B2 (en) * | 2003-08-28 | 2010-09-28 | Henry Ford Health System | Fuzzy logic guided inverse treatment planning |
US7496173B2 (en) * | 2003-12-12 | 2009-02-24 | University Of Western Ontario | Method and system for optimizing dose delivery of radiation |
US7529339B2 (en) * | 2003-12-12 | 2009-05-05 | University Of Western Ontario | Method and system for optimizing dose delivery of radiation |
US8422631B2 (en) * | 2007-12-07 | 2013-04-16 | Mitsubishi Heavy Industries, Ltd. | Radiation therapy planning apparatus and radiation therapy planning method |
US7945022B2 (en) * | 2008-03-12 | 2011-05-17 | Sun Nuclear Corp. | Radiation therapy plan dose perturbation system and method |
US8403821B2 (en) * | 2008-03-17 | 2013-03-26 | Mitsubishi Heavy Industries, Ltd. | Radiotherapy apparatus controller and radiation irradiating method |
US8325878B2 (en) * | 2008-05-08 | 2012-12-04 | The Johns Hopkins University | Real-time dose computation for radiation therapy using graphics processing unit acceleration of the convolution/superposition dose computation method |
US8682414B2 (en) * | 2009-11-27 | 2014-03-25 | Mitsubishi Heavy Industries, Ltd. | Radiation therapy apparatus control method and radiation therapy apparatus controller |
US20130102830A1 (en) * | 2010-06-22 | 2013-04-25 | Karl Otto | System and method for estimating and manipulating estimated radiation dose |
US20120041685A1 (en) | 2010-07-13 | 2012-02-16 | Vanderbilt University | System and method for estimating radiation dose and distribution using medium-dependent-correction based algorithms (mdc) |
US20150038762A1 (en) * | 2011-12-13 | 2015-02-05 | Gsi Helmholtzzentrum Fuer Schwerionenforschung Gmbh | Method for preparing an irradiation plan for a moving target volume without a movement compensation |
US20150199457A1 (en) * | 2014-01-14 | 2015-07-16 | Mitsubishi Electric Research Laboratories, Inc. | System and Method for Planning a Radiation Therapy Treatment |
US20150202464A1 (en) * | 2014-01-23 | 2015-07-23 | Mitsubis | Multi-Criteria Optimization in Particle Beam Dose Optimization |
Non-Patent Citations (1)
Title |
---|
Mihaylov et al. "Analytic IMRT dose calculations utilizing Monte Carlo to predict MLC fluence modulation," NIHPA, Med Phys. Apr. 2006: 33 (4): 828-839. |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US11147985B2 (en) * | 2017-02-02 | 2021-10-19 | Koninklijke Philips N.V. | Warm start initialization for external beam radiotherapy plan optimization |
US11541252B2 (en) * | 2020-06-23 | 2023-01-03 | Varian Medical Systems, Inc. | Defining dose rate for pencil beam scanning |
Also Published As
Publication number | Publication date |
---|---|
US20150199457A1 (en) | 2015-07-16 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
US11633626B2 (en) | Methods for real-time image guided radiation therapy | |
US9251302B2 (en) | System and method for planning a radiation therapy treatment | |
US8401148B2 (en) | Non-voxel-based broad-beam (NVBB) algorithm for intensity modulated radiation therapy dose calculation and plan optimization | |
US10369381B2 (en) | System and method for robust intensity-modulated proton therapy planning | |
Papagiannis et al. | Current state of the art brachytherapy treatment planning dosimetry algorithms | |
Bauer et al. | Integration and evaluation of automated Monte Carlo simulations in the clinical practice of scanned proton and carbon ion beam therapy | |
Inaniwa et al. | Implementation of a triple Gaussian beam model with subdivision and redefinition against density heterogeneities in treatment planning for scanned carbon-ion radiotherapy | |
Yepes et al. | A GPU implementation of a track-repeating algorithm for proton radiotherapy dose calculations | |
US9067064B2 (en) | Optimization process for volumetric modulated arc therapy | |
JP2014528797A (en) | Dose calculation using inhomogeneity-compensated superposition for radiation therapy | |
Lu et al. | Fluence-convolution broad-beam (FCBB) dose calculation | |
Fujimoto et al. | GPU-based fast pencil beam algorithm for proton therapy | |
Penfold | Image reconstruction and Monte Carlo simulations in the development of proton computed tomography for applications in proton radiation therapy | |
WO2012037472A2 (en) | Gpu-based fast dose calculator for cancer therapy | |
US9486644B2 (en) | Method and system for dose determination of radiation therapy | |
Vazquez et al. | A deep learning-based approach for statistical robustness evaluation in proton therapy treatment planning: a feasibility study | |
Cufflin et al. | An investigation of the accuracy of Monte Carlo portal dosimetry for verification of IMRT with extended fields | |
JP2013248133A (en) | Information processing device, information processing method, and program | |
Dickmann et al. | Fluence-modulated proton CT optimized with patient-specific dose and variance objectives for proton dose calculation | |
US11957934B2 (en) | Methods and systems using modeling of crystalline materials for spot placement for radiation therapy | |
Andersson | Mathematical optimization of radiation therapy goal fulfillment | |
Sundström | Scenario dose calculation for robust optimization in proton therapy treatment planning | |
Bedford | A discrete ordinates Boltzmann solver for application to inverse planning of photons and protons | |
Scholz | Development and evaluation of advanced dose calculations for modern radiation therapy techniques | |
Ganesh | Teaching Session |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
AS | Assignment |
Owner name: MITSUBISHI ELECTRIC RESEARCH LABORATORIES, INC., M Free format text: ASSIGNMENT OF ASSIGNORS INTEREST;ASSIGNOR:BRAND, MATTHEW;REEL/FRAME:033342/0750 Effective date: 20140708 |
|
STCF | Information on status: patent grant |
Free format text: PATENTED CASE |
|
MAFP | Maintenance fee payment |
Free format text: PAYMENT OF MAINTENANCE FEE, 4TH YEAR, LARGE ENTITY (ORIGINAL EVENT CODE: M1551); ENTITY STATUS OF PATENT OWNER: LARGE ENTITY Year of fee payment: 4 |
|
AS | Assignment |
Owner name: MITSUBISHI ELECTRIC CORPORATION, JAPAN Free format text: ASSIGNMENT OF ASSIGNORS INTEREST;ASSIGNOR:MITSUBISHI ELECTRIC RESEARCH LABORATORIES, INC.;REEL/FRAME:049241/0815 Effective date: 20190520 |
|
MAFP | Maintenance fee payment |
Free format text: PAYMENT OF MAINTENANCE FEE, 8TH YEAR, LARGE ENTITY (ORIGINAL EVENT CODE: M1552); ENTITY STATUS OF PATENT OWNER: LARGE ENTITY Year of fee payment: 8 |
|
AS | Assignment |
Owner name: HITACHI, LTD., JAPAN Free format text: NUNC PRO TUNC ASSIGNMENT;ASSIGNOR:MITSUBISHI ELECTRIC CORPORATION;REEL/FRAME:067122/0371 Effective date: 20240305 |
|
AS | Assignment |
Owner name: HITACHI HIGH-TECH CORPORATION, JAPAN Free format text: ASSIGNMENT OF ASSIGNORS INTEREST;ASSIGNOR:HITACHI, LTD.;REEL/FRAME:068292/0514 Effective date: 20240401 |