Optimization methods for radiation therapy planning

Abstract

An optimization technique for use with radiation therapy planning that combines stochastic optimization techniques such as Particle Swarm Optimization (PSO) with deterministic techniques to solve for optimal and reliable locations for delivery of radiation doses to a targeted tumor while minimizing the radiation dose experienced by the surrounding critical structures such as normal tissues and organs.

Claims

1. A method for providing a radiotherapy or radiosurgery treatment plan, comprising: using a particle swarm as a solution to an optimization problem, wherein the solution is a location of the particle swarm in a minimum energy state; and providing the location for a radiation dose from a radiation source.

2. The method of claim 1, wherein the optimization problem is: $\begin{array}{c}\mathrm{minimize}\\ t\end{array}.Math..Math.{\stackrel{.}{D}}_j{t}_j-D^{*}.Math.\mathrm{subject}\mathrm{to}{t}_{j}0,$ where D* is an ideal dose distribution or a prescribed dose, {dot over (D)}.sub.j is a dose rate contributed by a j-th candidate beam, and t.sub.j is a beam-on time for the j-th candidate beam, a constraint t.sub.j0 is a non-negative beam-on time.

3. The method of claim 2, further comprising the step of: solving the optimization problem when two or more beam-on times t.sub.j are found so that a created dose distribution .sub.j{dot over (D)}.sub.jt.sub.j is as close to the ideal dose distribution or the prescribed dose D* as possible.

4. The method of claim 1, further comprising the step of: determining a length of time the radiation dose is transmitted to the location using one or more algorithms selected from the group consisting of: least squares, non-negative least squares, least distance programming, and linear programming.

5. The method of claim 1, wherein the using step further comprises the steps of: modeling a tumor as a geometric volume comprising geometric objects with potential functions; modeling the radiation source as a set of particles subject to forces from the potential functions from the geometric objects; achieving a steady location of the set of particles; translating a position and a trajectory of each particle of the set of particles to a location of the radiation source; and providing a treatment plan comprising the location of the radiation source.

6. The method of claim 5, further comprising the step of: determining a length of time the radiation dose is transmitted to the location using one or more algorithms selected from the group consisting of: least squares, non-negative least squares, least distance programming, and linear programming.

7. The method of claim 5, wherein each particle p.sub.i of the set of particles includes a location: {right arrow over (x)}.sub.l for each particle p.sub.i=[x.sub.i, y.sub.i, z.sub.i], a velocity: {right arrow over (v)}.sub.l for each particle p.sub.i=[v.sub.x.sub.i, v.sub.y.sub.i, v.sub.z.sub.i], a type: static or kinetic, and a potential function.

8. The method of to claim 7, further comprising the step of: pre-computing a static potential field U.sub.s for each particle p.sub.i of the set of particles that is of the static type.

9. The method of claim 7, wherein a total external potential incident on each particle p.sub.i of the set of particles is:
{right arrow over ()}U.sub.i={right arrow over ()}(U.sub.K-{p.sub.i.sub.}+U.sub.S) wherein U.sub.S is a static potential field of the particle p.sub.i of the set of particles that is of the static type and U.sub.K is a kinetic potential field of the particle p.sub.i of the set of particles that is of the kinetic type.

10. The method of claim 7, further comprising the step of: calculating an acceleration {right arrow over (a.sub.l)} of each particle p.sub.i of the set of particles to update the velocity and the location of each particle p.sub.i of the set of particles in each time step t, repsectively:
{right arrow over (v.sub.l)}(t)={right arrow over (v.sub.l)}(t1)+{right arrow over (a.sub.1)}(t1).Math.t $\stackrel{->}{x_i}\left(t\right)=\stackrel{->}{x_i}\left(t-1\right)+\frac{\stackrel{->}{v_i}\left(t\right)+\stackrel{->}{v_i}\left(t-1\right)}{2}.Math.t.$

11. A method for Gamma Knife radiosurgery, comprising the steps of: assigning a potential function $\frac{1}{.Math.r}$ to each particle, wherein r is a distance from a particle to a voxel, and a is a constant scaling factor directly proportional to a dose kernel: initializing kinetic particles K at random positions with zero initial velocities or zero momentum; calculating potential functions associated to each particle using an approximate solution obtained from a constrained integer-linear problem: $\frac{\mathrm{minimize}}{y_i{Z}^{+}}.Math.a.Math.y-b.Math.$ $\mathrm{subject}\mathrm{to}{.Math.y.Math.}_{1}n$ wherein a is a row vector containing a tumor volume, each kernel can cover with a high-density dose, y is a column vector representing kernel spot sizes, b is a cardinality of a set of tumor voxels T, and n is a cardinality of the kinetic particles K; computing a potential field U.sub.s of static particles and a potential field U.sub.K of kinetic particles K; updating in each iteration, a location of each kinetic particle of the kinetic particles K based on forces from the potential field U.sub.K of kinetic particles K until the location of each kinetic particle converge; outputting a location of a spherical high dose volume for each kinetic particle.

12. The method of claim 11, further comprising the step of: calculating optimal dwell times using a non-negative least squares problem:
minimize .sub.j{dot over (D)}.sub.jt.sub.jD*.sub.2.sup.2 subject to t.sub.j0 where {dot over (D)}.sub.j is a dose rate contributed by a j-th candidate beam, t.sub.i is a dwell time for the j-th candidate beam of each spherical high dose volume, and D* is an optimal desired dose distribution.

13. The method of claim 11, wherein the kernel spot sizes are homogeneous.

14. The method of claim 13, further comprising the step of: creating a dose-volume histogram of homogeneous kernel spot sizes.

15. A method for high-dose rate (HDR) brachytherapy, comprising the steps of: providing a tumor modeled as a geometric volume; mapping each kinetic particle to a radiation source; imposing a set of initial conditions on each kinetic particle, the set of initial conditions comprising defined starting positions and initial velocities in a direction of each needle of a plurality of needles; recording a trajectory of each kinetic particle; and calculating a final dose delivered using a position of each kinetic particle as a center of the radiation source.

16. The method of claim 15, further comprising the steps of: modeling each needle as the trajectory of each kinetic particle with a potential function $\frac{1}{r^2},$ where r is a distance from a kinetic particle, and placing randomly static particles on a surface of the geometric volume of the tumor with a potential function $\frac{1}{r^2},$ where r is a distance from a static particle.

17. The method of claim 15, wherein the tumor is located on a prostate organ.

18. The method of claim 15, further comprising the step of: calculating a dwell time using a least-distance programming (LDP) problem:
Minimize x.sub.2.sup.2+.sub.maxs.sub.2.sup.2+.sub.mint.sub.2.sup.2
subject to Dxsb.sub.max
Dx+tb.sub.min wherein vector b.sub.min is a minimum dose for each voxel and vector b.sub.max is a maximum dose for each voxel, matrix D is a dose matrix of each discretized position along the plurality of needles, s and t are variable constraints to ensure feasibility, .sub.max and .sub.min are weighting variables.

19. The method of claim 15, further comprising the step of: determining a length of time the final dose is transmitted to a location using one or more algorithms selected from the group consisting of: least squares, non-negative least squares, least distance programming, and linear programming.

20. The method of claim 15, wherein the radiation source is an Ir-192 source.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The preferred embodiments of the invention will be described in conjunction with the appended drawings provided to illustrate and not to the limit the invention, where like designations denote like elements, and in which:

(2) FIG. 1 illustrates the pseudo-code according to one embodiment of the algorithm according to the invention.

(3) FIG. 2 is a schematic diagram illustrating the algorithm according to the invention.

(4) FIG. 3 is a flow chart of an exemplary optimization method for radiation therapy planning according to the invention.

(5) FIG. 4 illustrates a dose-volume histogram according to one embodiment of the invention.

(6) FIG. 5 illustrates a dose-volume histogram according to another embodiment of the invention.

(7) FIG. 6 illustrates a dose-volume histogram according to one embodiment of the invention.

(8) FIG. 7 illustrates a dose-volume histogram according to another embodiment of the invention.

DETAILED DESCRIPTION

(9) The invention is based on the idea that physical models can be used to guide an optimization algorithm consisting of the interactions of a system of stationary and moving particles, such that a desired behavior can be simulated in order to obtain an optimum maximum or minimum state. Specifically, the invention is inspired by particle swarm optimization that evolves in time until mapping its state of minimum energy to optimal locations of radiation sources in radiotherapy and radiosurgery. In particular, the invention deviates from the classical PSO algorithm where a particle represents a potential solution of the optimization function by using the whole swarm or particle distribution as the solution.

(10) Tumors, critical organs and other tissues are modeled as geometric volumes, whose surfaces have an associated potential function. The radiation source is modeled as kinetic particles subject to the forces from the potential functions from both the particles and the various geometric objects. The final configuration of the swarm of particles including their trajectories is the treatment plan. The term swarm is also referred to as map. Each particle includes a particular location or position such that the particles in their entirety create the swarm or map.

(11) The geometric characteristics of the treatment case are propagated throughout all the system so that a particle in the interior of the tumor can be aware of near-by critical organs or other tissues.

(12) The steps typically involved with radiation therapy treatment planning patient imaging, target definition/structure contouring, dose prescription, beam configuration optimization, plan generation, quality assurance assist in formally defining the characteristics of the algorithm according to the invention.

(13) Patient representation uses patient geometrical information described using set notation where V represents the set of all voxels in a CT or MRI scan of the patient, T represents the set of tumor voxels, C.sub.i represents the set of voxels for the i.sup.th critical structure. It can be noted that T and C.sub.i are pairwise disjoint. Finally, T and C.sub.i denote the voxels on the surface of the tumor and critical structures.

(14) The ideal dose distribution is a function that maps all voxels in V into a real interval prescribed by the physician, which corresponds to the minimum and maximum dose per voxel. The objective dose distribution is presented as:
D*:V.fwdarw.[a,b](Equation 1.0)

(15) For critical structures, the minimum dose is 0, while the maximum dose is determined by organ/tissue specific radiation tolerance. For the tumor, the minimum dose is determined by tumor control probability, and the gap between minimum and maximum dose, i.e., ba indicates the required dose uniformity.

(16) Dose calculation includes a particle, p.sub.i defined as an entity that has a location: {right arrow over (x)}.sub.l=[x.sub.i, y.sub.i, z.sub.i], a velocity: {right arrow over (v)}.sub.l=[v.sub.x.sub.i, v.sub.v.sub.i, v.sub.z.sub.i], a type: static or kinetic, and a potential function (otherwise referred to as a potential field).

(17) Optimization is performed in two phases: PSO and deterministic optimization. For PSO, particles are divided in two disjoint sets according to their type K or S, for kinetic and static particles, respectively. The particle potential functions are created from the corresponding radiation source dose distribution and the prescribed dose distributions. For static particles, the static potential field, U.sub.s is pre-computed. For a kinetic particle p.sub.i, the total external force incident on it is computed by the negative gradient of the total potential field of all particles:
.sub.ji{right arrow over (F)}.sub.J={right arrow over ()}.sub.j=iU.sub.j=m.Math.{right arrow over (a)}.sub.l(Equation 1.2)
where m stands for the mass of particle p.sub.i. The total external potential incident on a particle can also be expressed by:
{right arrow over ()}U.sub.i={right arrow over ()}(U.sub.K-{p.sub.i.sub.}+U.sub.s)(Equation 1.3)
Each particle's acceleration, {right arrow over (a)}.sub.i, is calculated in order to update the velocity and position of each particle in a given time step t, respectively:

(18) $\begin{array}{cc}\stackrel{->}{v_i}\left(t\right)=\stackrel{->}{v_i}\left(t-1\right)+\stackrel{->}{a_i}\left(t-1\right).Math.t& \left(\mathrm{Equation}1.4\right)\\ \stackrel{->}{x_i}\left(t\right)=\stackrel{->}{x_i}\left(t-1\right)+\frac{\stackrel{->}{v_i}\left(t\right)+\stackrel{->}{v_i}\left(t-1\right)}{2}.Math.t& \left(\mathrm{Equation}1.5\right)\end{array}$

(19) This process is repeated until a maximum time, T.sub.max, is reached or a tolerance error, , is achieved. Once a steady location of particles is achieved, particle positions are translated to radiation sources' positions in the treatment plan. In the deterministic optimization phase, numerical techniques such as non-negative least squares, least distance programming, linear programming, etc., are used to further refine the plan produced by PSO. The pseudo-code according to one embodiment of the algorithm is shown in FIG. 1.

(20) Given an initial kinetic particle swarm and a precomputed potential, each particle position gets updated according to the total external forces. This process is repeated until a stable configuration of particles is reached or a minimum error is achieved. The algorithm returns the final stable swarm, or final map of particle locations.

(21) Depending on the application, kinetic particles are initialized with an initial momentum. For the initial PSO phase, all kinetic particles are looped over one at a time. The current particle is excluded from the current total potential field (U). The negative gradient of the potential field (VU) is computed and the position of the current particle is updated using the motion equations 1.2, 1.4, and 1.5 above. After updating all particles, the total distance variation of the whole particle distribution is determined; if it is less than a given threshold, the process stops, otherwise the process repeats. In an alternative embodiment, the process stops when a number of iterations reaches a certain threshold value. In order to compute a final dose distribution, dwell times are determined and negligible dwell times are filtered out. A dwell time is the length of time for delivering a radiation dose. A dose-volume histogram can be generated from the final dose distribution. A dose-volume histogram facilitates understanding the quality of a treatment plan in terms of its dose distributions. Specifically, it provides the amount of volume per biological structure that received at least a certain percentage of the maximum dose.

(22) FIG. 2 is a schematic diagram illustrating the algorithm 100 according to the invention. A tumor mass 110, along with other critical structures 120 such as critical organs, and normal tissues are defined and modeled such as geometric volumes. The tumor mass 110 and other critical structures 120 are covered with static particles 130, each having its own potential function. It is contemplated that each potential function may be the same or different. More specifically, each potential function of each static particle 130 provide the potential field for the surface of the tumor volume and organ or tissue volumes.

(23) A radiation source is modeled as kinetic particles 140, or spherical high dose volumes, which are positioned within a cross-section of the tumor mass 110 or outside the tumor mass 110, but nearby. The kinetic particles 140 each have its own potential function. It is contemplated that each potential function may be the same or different. The kinetic particles 140 are subject to the forces from the potential functions from both the kinetic particles 140 and the various geometric objects. The kinetic particles 140 search for a minimum potential and stabilize.

(24) As shown in FIG. 2, the kinetic particles 140 are displaced toward the left side due to the placement of other critical structures 120 on the right side of the tumor mass 110. This effect is due to the potential field associated with the other critical structures 120. Specifically, the kinetic particles 140 move under the influence of a potential field. When the kinetic particles 140 are given an initial velocity (otherwise referred to as momentum) as shown by the arrow 150, the kinetic particles 140 traverse the tumor mass 110 in evenly spread trajectories 160. The final configuration of the swarm of particles, or spherical high dose volumes, including their trajectories 160 is the treatment plan.

(25) FIG. 3 illustrates a flow chart of the steps of an optimization method 200 for radiation therapy planning according to one embodiment of the invention. Structures are defined as one or more tumor masses 110 and surrounding critical structure 120 or normal structures. Static particles 130 are positioned on the surfaces of the critical structures 120 at step 202. The static particles 130 are positively charged creating a static potential field (as defined by a function). It is contemplated that the static potential field is constrained to a prescribed dose distribution. The static particles 130 may be applied randomly or according to a specified pattern.

(26) At step 204, one or more kinetic particles 140 or spherical high dose volumes are positioned near the tumor mass 110 or within the interior volume of the tumor mass 110, i.e., a cross-section of the tumor mass 110. The kinetic particles 140 are in an initial position such that all kinetic particles 140 define a map of locations, or swarm. An arbitrary potential field is applied to the kinetic particles 140 creating a dynamic potential field (as defined by a function). It is contemplated that the dynamic potential field is constrained to a dose distribution of a radiation source.

(27) At step 206, an initial velocity/momentum 150 is applied to each kinetic particle 140 of the swarm. The initial velocity/momentum 150 may be applied such that it is parallel to a principal implant direction of a radiation source. A force value incident upon the kinetic particles 140 is calculated from the potential fields. For example, the force may be a negative gradient of both the dynamic potential field and the static potential field. The force value incident upon each spherical high dose volume is used to arrive at an acceleration value.

(28) At step 208, the first location of each kinetic particle 140 is updated such that all kinetic particles 140 define a new map of locations. The kinetic particles 140 continue to move until the maps converge to a variation of distance that is less than a given threshold value at step 210. Alternately, the kinetic particles 140 move until a number of iterations reaches a certain threshold value.

(29) At step 212, the treatment plan is defined. The location of each kinetic particle 140 may be recorded or simply translated to define the radiosurgery plan. A radiation source executes the treatment plan. Specifically, location of each kinetic particle 140 translates to a radiation dose that is generated by a radiation source. Depending on the modality used, the point for receiving a radiation dose may be a beam location or interstitial implant trajectories 160.

(30) The optimization method 200 for radiation therapy planning according FIG. 3 may further comprise determining a length of time the radiation dose is to be delivered (i.e., dwell time) such that the dose distribution is within a prescribed ideal dose or at least as close to a prescribed ideal dose as possible. The length of time the radiation dose is to be transmitted may be determined using one or more algorithms selected from the following: least squares, non-negative least squares, least distance programming, linear programming, etc. The optimization method 200 according to the invention may also be used to create a dose-volume histogram from a dose distribution.

(31) According to one particular embodiment, the invention is applied to Gamma Knife radiosurgery. Dose kernels are mapped to a potential function associated to each kinetic particle 140 in the swarm and the interaction among them converges to a stable configuration. The locations of the kinetic particles 140 are then translated to dwell locations and the non-negative least squares algorithm is used to calculate the dwell time of each kernel aiming to meet the dose prescription.

(32) Each kinetic particle 140 is assigned a potential function

(33) $\frac{1}{.Math.r},$
where r is the distance from the kinetic particle 140 to a voxel, and a is a constant scaling factor directly proportional to the dose kernel radio to be used in the simulation. The intuition behind such an arrangement is that the kinetic particles 140, i.e. the spherical high dose volumes, should not be too dose to each other, and should not be too close to the boundary or the critical structures 120.

(34) The kinetic particles 140 are initialized at random positions with zero initial velocities or zero momentum 150. The quantity and potential functions associated to each kinetic particle 140 are calculated from an approximate solution obtained from the following constrained integer-linear problem:

(35) $\frac{\mathrm{minimize}}{y_i{Z}^{+}}.Math.a.Math.y-b.Math.$$\mathrm{subject}\mathrm{to}{.Math.y.Math.}_{1}n$
where a is a row vector containing the volume each kernel can cover with a high density dose, y is a column vector representing the distinct kernel spot sizes to be used, b is the cardinality of T, and n is the cardinality of K. The potentials U.sub.s and U.sub.K created by static and kinetic particles 140 respectively are computed. In each iteration, the locations of the kinetic particles 140 are updated based on the forces from the potential field until they converge. The PSO phase outputs the locations of each spherical high dose volume. To complete a particular radiosurgery plan, the dwell or beam-on times may be needed for each location. With {dot over (D)}.sub.j and t.sub.j be the dose rate and beam-on times for each particle/spherical high dose volume, and D* be the optimal desired dose distribution. Non-negative least squares is used to calculate the optimal beam-on times:
minimize .sub.j{dot over (D)}.sub.jt.sub.jD*.sub.2.sup.2 subject to t.sub.j0

(36) The optimization algorithm as applied to a C-shaped tumor 3D geometric object surrounding a spherical critical structure includes a prescription targeting 100 arbitrary dose units to the tumor while delivering 0 arbitrary dose units to the critical organ, which achieved convergence within 20 to 40 iterations.

(37) Homogeneous kernel spot sizes were used. Kinetic particles start at random locations inside the tumor and at each iteration they evolve, spreading evenly throughout the volume through the dynamic interactions among them and the potential field from the structure's surfaces. Particles get pushed in the opposite direction from the critical organ, which reflects the geometrical awareness imposed onto them by the potential field associated with the critical organ. Non-negative least squares is used to filter out those locations that are either redundant or that do not contribute to meeting the prescription. The particles distribute evenly along the tumor's volume while avoiding regions close to the critical organ.

(38) FIG. 4 and FIG. 5 illustrate dose-volume histograms. Specifically, FIG. 4 illustrates a dose-volume histogram of homogeneous kernel spot sizes using 4 mm dose kernels. FIG. 5 illustrates a dose-volume histogram of heterogeneous kernel spot sizes using 4 mm and 8 mm dose kernels. In the histograms, the horizontal axis is the dose, while the vertical axis is the percent volume. There is a line for each structure of interest, which describes the amount of dose delivered to a given percent volume. Specifically, line A is tissue structure, line B is tumor structure, line C is critical structure (i.e., organ), and line D corresponds to the prescription. The quality of the plan is comparable to the manually obtained clinical plans, which aims to cover the target tumor volume with 50% of the maximum dose.

(39) According to another particular embodiment, the invention is applied to high-dose rate (HDR) brachytherapy. As mentioned previously, HDR brachytherapy is performed as a four-step process: catheter or implant placement, image contouring, dwell time optimization, and dose delivery.

(40) Each kinetic particle is mapped to a radiation source and a set of initial conditions is imposed, which include defining starting positions and initial velocities in the direction of previous clinical implants. The key idea of this embodiment is to let the particles in the swarm to traverse the tumor and record its particles' trajectories. The trajectory of each particle is mapped as the implant trajectory, and the final dose delivered is calculated using each particle position as the center of a radiation source, further optimized using least-distance programming.

(41) Each needle is modeled as the trajectory of a kinetic particle with a potential function

(42) $\frac{1}{r^2},$
where r is the distance from the particle. Static particles are randomly placed on the surface of the prostate and various critical structures such as the rectum, urethra, and bladder. The potential function of each static particle is

(43) $\frac{1}{r^2}.$
The key to this embodiment is the initial velocity of the kinetic particles. Clinically, the bevel needles used have a curvature constraint and can bend up to a certain degree value. The initial velocity of the kinetic particles ensures that the curvature constraint will not be violated during the simulation. The total number of kinetic particlesthe number of needlesis specified by the user.

(44) To determine the initial position of the kinetic particles, the kinetic particles are randomly placed within a cross-section of the tumor with zero velocity and constrained to move within that cross-section. Once the particles stabilize, they are given a velocity vector parallel to the principal needle direction typically used in a clinical setting. After the trajectories converge, a third-degree polynomial regression is applied to smooth the final particles' trajectories, which are the final needle positions. To calculate the dwell time, the following formulation of a least-distance programming (LDP) problem is used:
Minimize x.sub.2.sup.2+.sub.maxs.sub.2.sup.2+.sub.mint.sub.2.sup.2
subject to Dxsb.sub.max
Dx+tb.sub.min
The vectors b.sub.min and b.sub.max specify the minimum and maximum dose for each voxel. The matrix D is the dose matrix of each discretized position along the needles. Ideally, all prescription constraints are met and Dxb.sub.min and Dxb.sub.max, while minimizing the total treatment time x. However, this is never feasible. To overcome this, slack variables s and t are added to the constraints to ensure feasibility. The weighting variables .sub.max and .sub.min (both are diagonal square matrices) reflect the importance of the various constraints. The goal of the objective function is to minimize the total treatment time while making sure that as few as possible voxels are violating the prescription constraints.

(45) Given a specified number of implants, the algorithm according to the invention finds the optimal trajectories that meet clinical goals. In addition, the algorithm reduces the number of implants. It is contemplated that implant procedures may be performed under the guidance of such pre-calculated trajectories. If the pre-calculated implant trajectories are not followed perfectly, the impact is negligible. Furthermore, it is shown that the algorithm obtains homogeneous dose plans.

(46) FIG. 6 shows the dose-volume histogram comparisons between the plan achieved by the algorithm of the invention (denoted by capital letters A, B, C, D) and the previously used clinical plan (denoted by lowercase letters a, b, c, d). Each letter combination shown in FIG. 6 illustrates a particular organ: A, a represents the rectum, B, b represents the bladder, C, c represents the urethra, D, d represents the prostate.

(47) With the dose-volume histogram illustrating the amount of dose delivered to a given percent volume, line D for the prostate in denotes that 98% of the prostate receives at least the prescribed dose. As shown in FIG. 6, all the prescription goals are achieved and is comparable to the previously used clinical plan. Although the algorithm of the invention produces less uniform results as shown by lines A, B, C, D than the clinical plan shown by lines a, b, c, d this may represent a superior result in that these dose inhomogeneities are mainly localized inside the target, which may allow higher cell damage and tumor control, while meeting the radiation tolerance constraints for all critical structures.

(48) The algorithm proves to produce high quality plans for a variety of setups varying in the number of needles used. A higher implant density is able to lower dose impact to the urethra, rectum and bladder. Finally, while the clinical plan for this case used 17 needles, the invention shows to be able to produce high quality plans using 10 needles corresponding to a 41% reduction.

(49) It is expected that implant placement procedures may be performed under the guidance of pre-calculated trajectories. In order to assess the applicability of our generated plans, we randomly perturb the trajectories to mimic the manual implant errors. However, least distance programming is able to compensate for these errors in perturbations.

(50) Even with the understanding that HDR brachytherapy plans inherently have inhomogeneous dose distributions, due to the nature of the radiation sources that deliver a very high dose to those tissues close to the implants, the versatility of the optimization technique according to the invention is tested by trying to obtain uniform dose distributions. This is done by varying the weighting variables .sub.max and .sub.min to reflect the importance of lower doses deposited inside the target. Higher dose homogeneity can be obtained, but with a slight change in the objective function of the least distance programming optimization, i.e.:
Minimize (x.sub.2.sup.2)+(1) (.sub.maxs.sub.2.sup.2+.sub.mint.sub.2.sup.2)
subject to Dxsb.sub.max
Dx+tb.sub.min where 01

(51) FIG. 7 displays the results. Similar to FIG. 6, each letter combination shown in FIG. 7 illustrates a particular organ: A, a represents the rectum, B, b represents the bladder, C, c represents the urethra, D, d represents the prostate. As shown in FIG. 7, the capital letters A, B, C, D denote the plan achieved by the modified objective function and lowercase letters a, b, c, d denote the previously used algorithm (capital letters A, B, C, D of FIG. 6).

(52) As shown, higher implant density improves homogeneity, but since the algorithm according to the invention is able to generate implant trajectories that bend intelligently along the prostate boundary, the critical organs rectum and bladder always remain underdosed (with respect to the clinical plan), while the urethra suffers a slight overdosage due to a better implant utilization inside the prostate.

(53) As described above, the invention was applied to Gamma Knife radiosurgery and High-Dose Rate Brachytherapy; however, any radiation therapy modality is contemplated including for example a particle therapy machine or a lower dose rate Brachytherapy seed.

(54) While the disclosure is susceptible to various modifications and alternative forms, specific exemplary embodiments of the invention have been shown by way of example in the drawings and have been described in detail. It should be understood, however, that there is no intent to limit the disclosure to the particular embodiments disclosed, but on the contrary, the intention is to cover all modifications, equivalents, and alternatives falling within the scope of the disclosure as defined by the appended claims.