SYSTEMS AND METHOD FOR DEVELOPING RADIATION DOSAGE CONTROL PLANS USING A PARETOFRONT (PARETO SURFACE)

Abstract

System, method, and computer program product to select a desired portion of a subject to receive a radiation dose, including: determining a plurality of Pareto points on a Pareto surface (PS); selecting a first reference point (p.sub.1) on a back side of an assumed Pareto surface (PS′) and first direction (q.sub.1) emerging therefrom towards PS′ from behind; selecting a first starting adjustment (x.sub.10) and iteratively developing forward a minimum criterion in steps until a final adjustment (x.sub.11) is reached that still is implementable in the radiation apparatus; stopping the forward development, thereby determining x.sub.11 represented by a final front point (y.sub.11) as a real Pareto point of PS′; and along q.sub.1, dismissing undetermined portions of the objective space in front of and behind y.sub.11 as not containing parts of the PS, and continuing with other remaining more determined portions that are assumed to each contain a part of PS′.

Claims

1. A computer program product the computer program product comprising computer readable program code to select a desired portion of a subject to receive a radiation dose by determining or developing a previously unknown multi-dimensional surface shaped by Pareto optimal points, situated in the objective space, the Pareto optimal points defining adjustments of a given radiation apparatus, when the code is executed on a computer, the executing code for . . . a) determining a plurality of Pareto points on a Pareto surface (PS) using a plurality of reference points and development directions, emerging therefrom; wherein b) selecting a first reference point (p.sub.1) on a back side of an assumed Pareto surface (PS′) and first direction (q.sub.1) emerging therefrom towards the assumed Pareto surface from behind; c) selecting a first starting adjustment (x.sub.10) and iteratively developing forward a minimum criterion in steps until a final adjustment (x.sub.11) is reached that still is implementable in the radiation apparatus(60,62,69); and stopping the forward development, thereby determining the achieved final adjustment (x.sub.11) represented by a final front point (y.sub.11) as a real Pareto point of the assumed Pareto surface (PS′); d) along the direction (q.sub.1), dismissing undetermined portions of the objective space in front of and behind of the achieved final front point (y.sub.11) as not containing parts of the PS, and continuing with other remaining more determined portions that are assumed to each contain a part of the assumed Pareto surface (PS′).

2. The computer program product of claim 1, the apparatus adjustment is optimized by y=F (x) and its image under the evaluation criteria is a point of the real Pareto surface (PS).

3. The computer program product of claim 1, comprising computer readable program code for repeating steps b) through d) with a next reference point (p.sub.2) and a next direction (q.sub.2) to a next final adjustment (x.sub.21), represented by a next final front point (y.sub.21) on the Pareto surface (PS′), the next final adjustment still being implementable in the radiation apparatus (60,62,69).

4. The computer program product of claim 3, comprising computer readable program code for repeating steps b) through d) with a further reference point (p.sub.3) and a further direction (q.sub.3) to a further final adjustment (x.sub.31) that still is implementable in the radiation apparatus (60,62,69), by stopping a further forward development, thus determining the further final adjustment (x.sub.31), represented by a further final front point (y.sub.31) as a Pareto point of the real Pareto surface (PS), dismissing undetermined portions of the objective space in front of and behind the final front point (y.sub.31) along the third direction (q.sub.3)

5. The computer program product of claim 2, continuing with other remaining undetermined portions with still further reference points (p.sub.i) and still further directions (q.sub.i), whereby the other remaining undetermined portions are assumed to each contain a part of the assumed Pareto surface (PS′) to build the real Pareto surface (PS).

6. The computer program product of claim 2, comprising computer readable program code for repeating until a stopping criterion is met (103), wherein each repetition comprising step b) through step d) with a further reference point (p.sub.i), a further direction (q.sub.i) and a further starting adjustment to obtain a further final adjustment (x.sub.i1) that still is implementable in the radiation apparatus (60,62,69), represented by a further front point (y.sub.i1); to yield a further final front point (y.sub.i1) by each repetition, and gain a plurality of further final front points; the stopping criterion (103) being associated with an error measure of a remaining, more determined portion, determined by a length, size or a volume thereof, to build up a determined Pareto front as multi-dimensional Pareto surface of given approximation and thus approximation quality.

7. The computer program product of claim 1 comprising computer readable program code on a tangible non-transitory medium.

8. The computer program product of claim 1 comprising code for determining reference points and directions including a reference point being the upper corner and the direction being a diagonal of an m-dimensional box containing part of the assumed Pareto Surface (PS′).

9. The computer program product of claim 1 comprising code for iteratively developing forward a minimum criterion in steps by solving a Pascoletti-Serafini problem along the direction that emerges from the starting point.

10. The computer program product of claim 1 wherein the kept portions are a union of m-dimensional boxes.

11. The computer program of claim 8, wherein the stopping criterion is met if, for some measure of a box, the obtained value for the largest box falls below a threshold value.

12. A method to determine or develop a previously unknown multi-dimensional surface shaped by Pareto optimal points, situated in the objective space, the Pareto optimal points defining adjustments of a radiation apparatus (60,62,69), the method comprising the steps a. determining Pareto points for a Pareto surface (PS); whereby b. selecting a first reference point (p.sub.1) on a back side of an assumed, not yet defined Pareto surface (PS′) and selecting a first direction (q.sub.1) emerging therefrom towards the assumed Pareto surface (PS′) from behind; c. selecting a first starting adjustment (x.sub.10) and iteratively developing forward along the first direction (q.sub.1) a minimum criterion until a final adjustment (x.sub.11) is reached that still is implementable in the radiation apparatus (60,62,69); and stopping the forward development, thereby determining the achieved final adjustment (x.sub.11) represented by a final front point (y.sub.11) as a real Pareto point of the assumed Pareto surface (PS′); d. dismissing undetermined portions of the objective space in front of and behind of the achieved final front point (y.sub.11) as not containing parts of the real Pareto surface.

13. The method of claim 12, the apparatus adjustment (70) is optimized by y=F (x) and its image under the evaluation criterion is a point of the real Pareto surface (PS).

14. The method of claim 12, repeating steps b) through d) with a next reference point (p.sub.2) and a next direction (q.sub.2) to a next final adjustment (x.sub.21), represented by a next final front point (y.sub.21) on the Pareto surface (PS′), the next final adjustment still being implementable in the radiation apparatus (60,62,69).

15. The method of claim 14, comprising computer readable program code for repeating steps b) through d) with a further reference point (p.sub.3) and a further direction (q.sub.3) to a further final adjustment (x.sub.31) that still is implementable in the radiation apparatus, and stopping the further forward development, thus determining the further final adjustment as next final adjustment (x.sub.31), represented by a further final front point (y.sub.31) as a Pareto point of the real Pareto surface (PS), dismissing undetermined portions of the objective space in front of and behind the further final front point (y.sub.31) along the third direction (q.sub.3).

16. The method of claim 15, wherein in step (c) of claim 12 the minimum criterion is repeatedly developed forward by solving Pascoletti-Serafini problems.

17. The method of claim 12, wherein in step (c) the minimum criterion is developed forward by solving a Pascoletti-Serafini problem.

18. The method of claim 12, repeating steps (b) to (d) until a stopping criterion is met, each repetition comprising step b) through step d) with a further reference point (p.sub.i), a further direction (q.sub.i) and a further starting adjustment to obtain a further final adjustment (x.sub.i1) that still is implementable in the radiation apparatus, represented by a further front point (y.sub.i1); to yield a further final front point (y.sub.i1) by each repetition, and gain a plurality of further final front points; the stopping criterion being associated with an error measure of a remaining, more determined portion, determined by a length, size or a volume thereof, to build up a determined Pareto front as multi-dimensional Pareto surface of given approximation and thus approximation quality.

19. The method of claim 13, wherein at least a few calculations along a few directions are associated with separate kernels of a multi kernel processor (80) to provide parallel processing.

20. The method of claim 12, wherein a plan is Pareto optimal if there is no other plan with at least equally low values for all objectives and a lower value for at least one objective and an image set of the Pareto optimal plans as points is the Pareto surface (PS). a. selecting a number of surface points of a Pareto surface (PS) that define m objective functions, and an initial m-dimensional starting box (B0) containing the surface points; b. computing a set of outer knee points and a set of inner knee points from the initial starting box (B0) and the surface points ( . . . ) to determine a set of m-dimensional boxes (B1 to B4) whose union covers at least a portion of the Pareto surface; c. increasing the number of surface points of the Pareto surface by adding a point in one of the boxes (B1 to B4) whose union covers the portion of the Pareto surface.

Description

INTRODUCTION TO THE DRAWINGS

[0057] Various other features of the present invention will be made apparent from the following detailed description and the drawings relating to embodiments. They are seen to be exemplary only, not being read as limitations into the claims.

[0058] FIG. 1 shows an example of a workflow of the disclosed method.

[0059] FIG. 1a shows a radiation apparatus 60, 62 and 69, taken from U.S. Pat. No. 7,391,026.

[0060] FIG. 2 exemplifies how radiation dose affects multiple volumes of interest (VOIs). With each VOI a dose evaluating objective function ƒ.sub.i is associated, measuring some dose average.

[0061] FIGS. 3a to 3d explain how iteratively a better approximation of a Pareto surface PS' is attained, e.g. by cutting away parts that cannot contain the Pareto surface PS. In this way a covering of the (real) Pareto surface (PS) is obtained by “boxes” (2- or 3- or multi-dimensional). The cut away boxes are 51a, 51b in FIGS. 3b; 5, 5 and 52a, 52b in FIG. 3d.

[0062] FIGS. 3e to 3h further explain FIG. 3d with respect to inner knee points and outer knee points as well as their sets and the sets of boxes B1, B2, B3 and B4 that knee points define. The initial box used for starting is shown stippled in FIG. 3h (and FIGS. 3f and 3g).

[0063] FIG. 4 shows multiple reference points and directions that can be solved simultaneously. Thereby obtained points on the Pareto surface PS enable a construction of a covering set of boxes.

[0064] FIGS. 5a and 5b display Pareto surface approximations, obtained by hyperboxing . . .

[0065] FIG. 5a shows the points with the ideal point (0,0,0) in front

[0066] FIG. 5b shows the triangulated points seen from a side.

[0067] FIGS. 6a to 6e are alternative Pareto surface examples of triangulated Pareto surfaces PS as developed by examples of the disclosed method. FIGS. 6a to 6e show the PS from different sides.

[0068] FIG. 7 shows a computer with multiple cores or kernels used in the methods disclosed herein.

DETAILED DESCRIPTION OF THE EMBODIMENTS

[0069] One approach is the planning problem in intensity modulated radiation therapy (IMRT) and in volumetric arc therapy (VMAT).

[0070] The radiation apparatus is not shown, as suitable ones can be found in the prior art for IMRT and VMAT therapy.

[0071] In these planning processes, a compromise between target coverage (high radiation does) and sparing (low radiation dose) of nearby critical volumes (or structures) will be found. Therefore, this constitutes an optimization problem with multiple conflicting objectives (MOP). A MOP can be represented as

min.sub.xεX{ƒ.sub.i(x)},=1, . . . ,m

where

[0072] x is the decision variable

[0073] X denotes the feasible set

[0074] ƒ.sub.i(x) are the objectives such that smaller values denote preferable properties of x.

[0075] In the radiation therapy context, the decision variables are the plan parameters. In particular, x may denote the fluence profile of the beam or arc fields. A description of the MOP for IMRT can be found in Kuefer et al. and is incorporated herein by reference as to the definition of a Pareto multi-dimensional surface named “PS”.

[0076] A Pareto point is a point carrying the definition of a radiation plan and being of such nature that no quality indication on the axes of the point can be improved without deteriorating a quality indication on any other of the given axes.

[0077] A plan is Pareto optimal if there is no other plan with at least equally low values for all objectives and a strictly lower value for at least one objective. The image set (or group) of the Pareto optimal points (POPS) is referred to as the Pareto surface (PS).

[0078] The examples target this MOP by computing, using a database holding Pareto optimal plans which may then be “explored” (used for evaluation of a suitable plan).

[0079] In that the invention allows the exploration of all available compromises, it differs from a range of other multi-objective optimization strategies which require specifying of the planner's preferences beforehand, such as solving weighted sum problems with fixed weights, goal programming, and lexicographic ordering.

[0080] In that the invention allows the exploration of non-convex parts of the Pareto surface, it differs from any weight-based scheme which, once or iteratively, solves weighted sum problems.

[0081] As stated above, the invention computes a database of plans on the Pareto surface (PS). This is a consequence of the fact that the PS cannot, in general, be determined in closed form. Rather, the PS must be approximated by individual points.

[0082] Based on the points computed in previous iterations, an upper and lower bound on the region where the PS is situated is provided, allowing for targeted calculation of points in places where the Pareto surface PS is not yet represented well enough.

[0083] From the representation points, a surface approximation of the Pareto surface PS can be obtained by triangulation schemes, such as the Delaunay triangulation.

[0084] Individual solutions to the multi-objective problem can be computed from scalarizations.

[0085] A most common scalarization method is the weighted sum scalarization method. However, this method does not compute solutions in non-convex parts of the PS. Thus it is not suitable to calculate a representation set in a non-convex setting, as arbitrarily large areas or parts of the Pareto surface PS are missed.

[0086] The Pascoletti-Serafini method as introduced earlier overcomes this drawback. By varying its parameters, it can obtain all Pareto optimal points on an arbitrary PS, in particular in convex and non-convex parts of a given PS alike. The Pascoletti-Serafini scalarization solves problems of the form

α.fwdarw.min

s.t. p+αq=F(x)+λ

(α,x,λ)εR×X×R.sub.+.sup.m

Where . . .

F(x)=(ƒ.sub.1(x),ƒ.sub.2(x), . . . ,ƒ.sub.m(x))

is the vector of objective function values and p and q are the parameters.

[0087] For each Pareto optimal point x there is a parameter set (p,q) such that x is a solution to the corresponding Pascoletti-Serafini problem. We refer to p as the reference point and to q as the direction of the Pascoletti-Serafini problem, see FIG. 3 explaining the naming.

[0088] The fact that a non-convex Pareto surface PS or non-convex portion of a Pareto surface PS can be represented with solutions from differently parametrized Pascoletti-Serafini problems is used in the context of radiation therapy planning. There, plan quality is frequently assessed by dose and volume quantiles such as the target volume percentage receiving the prescribed dose, or the dose achieved in the target's 95% quantile. Translating these clinical goals into objectives yields non-convex functions and a non-convex PS.

[0089] The disclosure describes an iterative determination of the parameters p and q for the Pascoletti-Serafini problem such that an economic representation of an arbitrary Pareto surface is obtained. M-dimensional boxes are found and used, such that the Pareto surface PS is contained in their union. The box dimensions provide a measure for the approximation quality realized by the previously calculated points. The box dimensions also provide the means to determine where a new point on the Pareto surface PS should be calculated, and what is used as best parameters p and q to obtain that point.

[0090] Because of this usage of m-dimensional boxes, we refer to the algorithm as “hyperboxing”.

[0091] The idea and concept of hyperboxing starts with an m-dimensional box containing the Pareto surface PS. In a first iteration the method solves a Pascoletti-Serafini scalarization problem “along a diagonal of the box” Thereby a minimal vertex of the box is used as the reference point p and the diagonal of the box from minimal to maximal vertex is used as direction q, see FIG. 3a.

[0092] By solving the Pascoletti-Serafini problem with the parameters p and q, a solution is found which is added to the representation set. Then the initial box is removed from the set of covering boxes, and new boxes are added to the set of covering boxes. These have been constructed from the initial box and the representation point found before.

[0093] In a following iteration the largest box from the covering boxes is chosen (by any measure, used to identify the size of a box). Within this box, the Pareto surface PS is represented most poorly so far. The Pascoletti-Serafini scalarization along the diagonal of the chosen box is solved, the solution is added as a point to the representation set, and the set of covering boxes is updated again. This iteration is repeated, until a stopping criterion (or stop criterion) with respect to a size of the largest box is met.

[0094] Applied to FIG. 3b, the following iteration chooses the largest box 51a from the covering boxes (by any measure, used to identify the size of a box). Within this box, the Pareto surface PS is represented most poorly so far. The Pascoletti-Serafini scalarization along the diagonal q.sub.1 of the chosen box is solved, the solution is added as a point y.sub.11 to the representation set, and the set of covering boxes is updated again. This applies in the next iteration with y.sub.31 and y.sub.21 in FIG. 3c, where the diagonals are q.sub.2 and q.sub.3. The iteration using the cut away boxes 5 and 52a is shown in FIG. 3d. Points y.sub.31 and y.sub.21 are added to the representation set. The set of covering boxes is then updated again.

[0095] This iteration is repeated, until a stopping criterion (or stop criterion) with respect to a size of the largest box is met.

[0096] The hyperboxing method maintains a set of boxes that, as a union, cover the Pareto surface PS. In each iteration step, larger boxes are replaced by smaller ones, strictly reducing the covered area and thus approximating the Pareto surface PS ever more tightly. This tightening is illustrated in FIG. 3, i.e., FIG. 3a to FIG. 3d, as mentioned before.

[0097] Assume having obtained (α, x.sub.11, λ) as the solution of the Pascoletti-Serafini problem along the diagonal q.sub.1 of an initial box, as depicted in FIG. 3a. This defines the Pareto point y.sub.11=F (x.sub.11). It holds that there are no Pareto optimal points in the area left-hand and below the point y.sub.11 (the point marked square, not round on the diagonal q in FIG. 3b). Otherwise a lower value for a could have been attained in the Pascoletti-Serafini problem. The apparatus adjustment is optimized by y=F (x) and its image under the evaluation criteria is a point of the real Pareto surface (PS).

[0098] Also, there are no Pareto optimal points in right-hand and above y.sub.11 (the point marked as square on the Pareto surface PS) by definition of Pareto optimality, as shown before. Thus, the Pareto surface PS is found to be contained in the union of two remaining boxes obtained by cutting away the hatched area from the initial box, see FIG. 3b.

[0099] The same method step is used in the next and further iterations, leading to a set of more and smaller boxes, approximating the Pareto surface PS more tightly as shown in FIG. 3c and FIG. 3d.

[0100] Several approaches to approximate a non-convex PS have been proposed in the literature. Among those are specific algorithms to approximate the non-dominated set of a non-convex bi-criterion problem. These have never been generalized to higher dimensions, as is needed in the context of radiation treatment planning. Such generalizations to higher dimensions (as needed in the context of radiation treatment planning) are non-trivial.

[0101] Among others, generating the non-inferior set in mixed integer bi-objective linear programs: An application to a location problem, see Solanki, page 1. He considers a rectangle representation, by defining an initial rectangle containing the non-dominated set and a finite set of non-dominated points, the rectangle representation comprises the set of rectangles spanned by neighbouring points such that the non-dominated set is contained in the union of the rectangles. Each rectangle is associated with an error, namely the maximum of the rectangle sides relative to the side lengths of the initial rectangle. This prior art algorithm thus successively solves an augmented weighted Tchebycheff-Problem along the diagonal of the box with the largest error.

[0102] For higher dimensions than two, un-biased approximation in a multicriteria optimization is suggested, by Klamroth et al., page 1. This describes both, an inner and an outer approximation algorithm for the approximation non-convex non-dominated sets. The principle of splitting boxes with a weighted Tchebycheff-Method is applied here: for each of the (inner respectively outer) approximation points, an opposite point is constructed, and the lexicographic weighted Tchebycheff-Problem is solved with weights corresponding to the box dimensions. The issue there is that determining the region where the approximation is worst is a difficult disjunctive problem. As a consequence, in each iteration of the algorithm the Tchebycheff-Problems are solved for all current boxes, not for one at a time as suggested by the methods of the disclosure, called “hyperboxing”.

[0103] In prior art this leads to long calculation times. The disclosure alleviates this drawback as it allows targeted refinement of the approximation on one box at a time. And using several Kernels of processors, each Kernel may calculate one box at a time, thus many boxes may still be calculated at a time, but distributed to more than one Kernel, see FIG. 7.

[0104] In contrast to the disclosure supplied herein, several other approaches of prior art do not systematically update an inner and outer approximation of the Pareto surface PS but choose each reference point and each direction a priori or randomly.

[0105] FIG. 1 shows one preferred workflow of hyperboxing. FIG. 1a shows a radiation apparatus as taken e.g. from prior art U.S. Pat. No. 7,391,026 (Fraunhofer) and its FIG. 1.

[0106] First, an initial starting box is defined at step 100. This box can be chosen such that only interesting regions of the Pareto surface PS are approximated. It can be chosen to be centred at a specifically chosen starting solution, or it can be large enough to contain the whole PS.

[0107] It may be beneficial to choose the initial starting box as small as possible, if the whole Pareto surface PS is to be contained. This will be achieved by e.g. computing the individual minima with respect to all objectives and choosing a smallest box that contains all these points as a starting box.

[0108] While it is possible that the Pareto surface PS exceeds this region, it is a reasonable estimate to start with.

[0109] Finding the true minimal starting box that contains the PS is a computational problem in dimensions greater than two.

[0110] In a next step 102, to be repeated iteratively, the “largest” box is found from the set of covering boxes. As measure for the size of a box, different metrics are possible. One possibility is to use the volume. A computationally more advantageous metric is the smallest side length of a box . . .

μ(B)=min{L.sub.i,i=1, . . . ,m}.

[0111] Then it is checked at 103 if the measure of the largest box went below a threshold value, acting as a stop(ping) criterion. If this is the case, the method finishes at 106. Otherwise, it proceeds to step 104 with searching for an additional point on the Pareto surface PS.

[0112] To find the point, the Pascoletti-Serafini-Problem along q.sub.i of the chosen box is solved. “i” is the iteration index of loop of steps 102 to 105.

[0113] In a next step 105, the set of covering boxes may be updated. This includes removing the previously largest box, and calculating the new boxes from the old boxes and the newly found solution, as shown in FIG. 3c.

[0114] This may, for example, be achieved by updating both, the set of so-called “outer knee points” and the set of so-called “inner knee points”. The outer knee points form a set of points with the following property: for each point on the Pareto surface PS there is an outer knee point that is smaller or equal in every component. Analogously, the inner knee points form a set with the following property: for each point on the Pareto surface PS there is an inner knee point that is greater or equal in every component, see J. Legriel, C. LeGuernic, S. Cotton, O. Mahler, “Approximating the Pareto Front of Multi-Criteria Optimization Problems”, Tools and Algorithms for the Construction and Analysis of Systems, Lecture Notes in Computer Science, Vol. 6015/2010, pages 69 to 83, 2010 (herein Legriel et al.).

[0115] The updated set of boxes is then defined by the set of pairs comprising one inner and one outer knee point from the updated sets, respectively, such that the two points span a box with non-zero volume.

[0116] The method then proceeds to 102 again with determining the largest box from the updated set of boxes at 105.

[0117] While it is possible to alternately find a new point and update the set of covering boxes, it is also possible to simultaneously solve multiple Pascoletti-Serafini problems at 104 before updating the set of boxes again.

[0118] This may be done by finding new points along the diagonals of the K≧2 largest boxes within each iteration, or by varying p and/or q in some way such that multiple points are searched for simultaneously, see FIG. 4 for an example.

[0119] This allows parallelization on a computer 80 with multiple cores (or multiple kernels) as shown in FIG. 7 and FIG. 1a. At least a few calculations along a few directions are associated with separate kernels of a multi-kernel processor 80. This multi-kernel processor 80 is depicted in FIG. 7. FIG. 7 is related to FIG. 4 that has three directions q.sub.1, q.sub.2, q.sub.3. The multi-kernel (multi core) processor provides parallel processing. The parallel processing is exemplified in FIG. 7, where a bus interface 85 is present that has a further bus 87, coupling the bus interface to the outside of the processor. The bus interface has bidirectional connection to four kernels (cores). Each kernel (core 1, core 2, core 3 and core 4) has a register file and an ALU. This is 81 for core 1, 82 for core 2 and 83 for core 3, finally displayed is 84 for core 4. Each register file has a bidirectional link to the ALU (in and out). The register file has a bidirectional connection to the bus interface, for each of the cores 1 to 4.

[0120] Applying FIG. 4 to FIG. 7 (or vice versa), calculation along q.sub.1 could be performed in core 1, calculation along q.sub.2 could be performed in core 2 and calculation for direction q.sub.3 could be performed in core 3. Thus at least a few calculations along a few directions q.sub.1, q.sub.2 and q.sub.3 are associated with separate kernels of the shown multi-kernel processor 80 to provide parallel processing.

[0121] Hence, the above-described provides for calculating well-placed points on the Pareto surface PS, and applying it to the IMRT or VMAT radiation planning problem.

[0122] The method may be implemented in software that is run by a computerized system 80. A feasible set and objective functions are used and this condition is satisfied by staying within a linear programming environment or within the context of a convex feasible set and convex objective functions. Several goals are obtained, that may be combined or used individually.

[0123] The computerized system 80 is schematically shown in FIG. 1a that also shows the radiation apparatus 60, 69 as taken from prior art U.S. Pat. No. 7,391,026, a patent that is also owned by the owner of this application. In summary, FIG. 1a shows on the right hand side the computerised system in a housing 90, having the computer 80 that may be built as FIG. 7 shows. The computer 80 has a bidirectional bus 87 that connects to an interface 70 that gives control signals x to the devices 60 and 69, which are part of the full radiation apparatus. The x-parameters are plans (comprising control variables that make up each plan) and the plan controls movable elements 62, 64 and 69 of the whole radiation apparatus.

[0124] Both device parts 60 and 69 have directions to move or operate. An axis z′ on the device 60 allows the turning a of device 62 with respect to a stationary part 60. The radiation beam head 64 can thus be turned along a.

[0125] The other part 69 of the radiation apparatus is not providing radiation beams, it is positioning the patient P. It can be positioned upwards and downwards along h and forward and backward along z (parallel to z′). This device may also be turned along an angle β (+80° to −40° are displayed). It receives along control line 71 signals, that make up the plan for providing the patient P with radiation doses that are as well determined by control signals x along control line 72, both emerging from interface 70. This interface receives its commands through bus 87, the commands emerge from the computerised system 80 using the PS surface.

[0126] Included in the housing 90 may also be a non-specifically displayed data base that is as well shown in the prior art radiation device of mentioned US '026 patent as database 1. The radiation apparatus 60, 62 and 69 is also explained in the corresponding European patent EP 1 438 692 B1, page 8, para [047] to [049]. Prior art radiation apparatus carries different reference signs, but has been included into this disclosure with other reference signs that more suitably match with the other reference signs used herein.

[0127] The computerized system 80 has a certain reading device to receive a computer program product, not displayed, either by inserting a readable device or system 80 is programmed by transferring the computer program product into the main memory (of the computerized system 80). The code that is transferred for the program to run is provided in this way into the computerized system 80.

[0128] When selecting a desired portion of a subject to receive a radiation dose provided by the beam head 64 in FIG. 1a, FIG. 2 shows this abstraction of four volumes of interest 1 to 4 as can be present in the patient P. Each volume has a dose evaluating objective function ƒ.sub.i, or is associated with it, measuring a dose average for the respective volume of interest. Thus the objective function ƒ.sub.1 is related to a dose average ƒ.sub.1. The same applies for volume of interest 2 which is associated with the objective function ƒ.sub.2, measuring the dose average for volume of interest 2. The same applies for the other two abstract volumes 3 and 4 with objective functions ƒ.sub.3 and ƒ.sub.4.

[0129] One of these volumes of interests may be the target, the other three ones may be the risks. Their location and distribution and size is taken from the whole body of the patient P who is to receive a treatment plan of radiation.

[0130] The hyperboxing covers the Pareto surface PS with m-dimensional boxes and iteratively approximates the Pareto surface to an arbitrary degree of exactness. The exactness may be specified by a certain size parameter that is larger or smaller, like a resolution error Δ when digitizing an analogue signal.

[0131] It may be specified, how to add a point to the Pareto surface PS to best improve the representation of the PS.

[0132] These benefits are not available in traditional approximation methods, such as the normal boundary intersection method, the epsilon constraint method, or the uniform weights scalarization approach. The benefits are not yet available in prior art for an at least partly or entirely non-convex Pareto surface PS, caused by non-convex objectives functions or non-convex feasible set.

[0133] For practical and computational reasons, the application may be limited to dimensions 6 (lower than seven dimensions). The computation of the inner and outer knee points and the implicit definition of the boxes as combinations of an inner and an outer knee point may limit dimensionality in practice, as well as the computational effort to calculate a Delaunay triangulation from the Pareto points found.

[0134] In step 105 of FIG. 2, the set of covering boxes may be updated (explained earlier). This includes removing the previously largest box, and calculating the new boxes from the old boxes and the newly found solution, as shown in FIG. 3c. The same method step is used in the next and further iterations, leading to a set of more and smaller boxes, approximating the Pareto surface PS more tightly as shown in FIG. 3c and FIG. 3d.

[0135] This may, for example, be achieved by updating both, the set of so-called “outer knee points” v.sub.i and the set of so-called “inner knee points” w.sub.i as present in FIG. 3d and highlighted in FIGS. 3f and 3g, both together in FIGS. 3e and 3h.

[0136] The outer knee points v of FIG. 3e form a set of points v.sub.1, v.sub.2, v.sub.3, v.sub.4, . . . with the following property: for each point (e. g. y.sub.11, y.sub.21, y.sub.31) on the Pareto surface PS there is an outer knee point that is smaller or equal in every component. Analogously, the inner knee points w of FIG. 3e form a set with the following property: for each point (e. g. y.sub.11, y.sub.21, y.sub.31) on the Pareto surface PS there is an inner knee point that is greater or equal in every component, see Legriel et al. earlier in this disclosure.

[0137] The updated set of boxes is then defined by the set of pairs comprising one inner knee point w.sub.i and one outer knee point v.sub.j from the updated sets, respectively, such that the two points span a box with non-zero volume.

[0138] When viewing at FIG. 3e, the outer knee points v and the inner knee points w of FIG. 3d can be easily spotted. In FIG. 3e they have names, v.sub.j (j is running from 1 to 4) and w.sub.i (i is running from 1 to 4). In FIG. 3d the names of the knee points are not yet presented. In both FIG. 3e and FIG. 3h the boxes are each spanned by a pair of knee points. Namely, these are box B1, spanned by knee points v.sub.1 and w.sub.2 (one outer and one inner knee point), box B2 is spanned by knee points v.sub.2 and w.sub.2, box B3 is spanned by knee points v.sub.3 and w.sub.3 and box B4 is spanned by knee points v.sub.4 and w.sub.4. These are boxes of non-zero volume that in their entirety form a set of covering boxes that approximate the Pareto surface PS.

[0139] Now separating FIG. 3e (or FIG. 3d) will yield FIGS. 3f and 3g. FIG. 3g shows the inner knee points w.sub.1, w.sub.2, w.sub.3, w.sub.4 and their relation to a respective point on the Pareto surface. FIG. 3f shows the outer knee points v.sub.1, v.sub.2, v.sub.3, v.sub.4, also related to a respective point on the Pareto surface. The points on the Pareto surface are

[0140] Using FIGS. 3g and 3f together and putting them on top of each other, will yield FIG. 3h, where the same structure is given as shown in FIG. 3d, outlining the boxes covering the Pareto front PS, as discussed before. This is exactly what is shown in FIG. 3d, just outlining or enhancing the points that are present in FIG. 3d already.

[0141] From FIGS. 3d, 3e it may be taken that from the initial box B0 (given by broken line in FIG. 3h) and the set of Pareto points y.sub.11, y.sub.21 and y.sub.31 a set of outer knee points v.sub.1, v.sub.2, v.sub.3, v.sub.4 and a set of inner knee points w.sub.1, w.sub.2, w.sub.3, w.sub.4 may be computed.

[0142] The initial box B0 is shown in broken lines in FIG. 3h (encompassing the set of boxes B1 to B4, each made up with a pair of knee points (outer and inner knee point). Thus a set of m-dimensional boxes may be determined. The union of these boxes covers the Pareto surface (at least a portion thereof). The number of points on the Pareto surface can now be increased by adding one point in one of the B1 to B4 boxes whose union covers the Pareto surface.

[0143] Thus surface points on the Pareto surface PS define m objective functions, using the initial m-dimensional starting box containing the points.

[0144] Disclosed is a system and a method for determining a desired portion of a subject to receive a radiation dose, including iteratively updating eligible reference points and directions to gradually build up a Pareto surface. By examining the previously calculated points on the Pareto surface and the eligible reference points and directions, a next reference point and direction is chosen, to yield a new Pareto optimal point (PoP) on the Pareto surface. This iteration step is repeated until a criterion to stop is met.

[0145] The invention has been described in terms of the various aspects and features, and it should be appreciated that many equivalents, alternatives, variations, and modifications, aside from those expressly disclosed, are in the skilled man's reach (PHOSITO's reach) and thus within the scope of the invention as claimed. Therefore, the disclosure should not be limited to a particular embodiment.