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 1 ) on a back side of an assumed Pareto surface (PS′) and first direction (q 1 ) emerging therefrom towards PS′ from behind; selecting a first starting adjustment (x 10 ) and iteratively developing forward a minimum criterion in steps until a final adjustment (x 11 ) is reached that still is implementable in the radiation apparatus; stopping the forward development, thereby determining x 11  represented by a final front point (y 11 ) as a real Pareto point of PS′; and along q 1 , dismissing undetermined portions of the objective space in front of and behind y 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′.

Description:
FIELD OF THE INVENTION 
       [0001]    This invention relates generally to program products, systems and methods for developing control plans for radiation dosage and, more particularly, to systems and methods to determine a desired plan for applying a radiation dose. 
       BACKGROUND OF THE INVENTION 
       [0002]    The application of radiation to an individual is commonplace in a number of industries, such as the medical industry, where radiation is used in imaging applications and therapy applications, and the security industry, where radiation is used to perform inspection processes. Regardless of the industry, whenever radiation is applied to a living subject, a balance is sought between the ability to achieve the best results of the process utilizing the radiation, which typically leans toward increased radiation dosage, and the safety of the subject, which typically leans toward decreased radiation dosage. 
         [0003]    The following references provide guidance: Intensity Modulated Radiotherapy: A Large Scale Multi Criteria Programming Problem, by authors Kuefer, Monz, Scherrer, Trinkaus, Bortfeld, Thieke, OR Spektrum, 2002; Unbiased approximation in multicriteria optimization (herein Kuefer et al.); Mathematical Methods in Operational Research, by Klamroth, Tind, and Wiecek, Vol. 56 (3), 2002, pages 413 to 437 (herein Klamroth et al.) and Approximation Methods in Multiobjective Programming. Journal of Optimization Theory and Applications, by Ruzika and Wiecek, Vol. 126, 2005, pages 473 to 501 (herein Ruzika et al.). Generating the non-inferior set in mixed integer bi-objective linear programs: An application to a location problem, by Solanki, European Journal of Operational Research, Vol. 68, 1993, pages 356 to 373 (herein Solanki). 
         [0004]    U.S. Pat. No. 8,489,366 B2 (Bortfeld, Grant) discloses a method for discovering a Pareto surface but will fail when the surface has at least one concave portion (bent inwardly). It operates according to triangles that are placed around a small segment of a Pareto surface, and then iterates further steps to further define the unknown surface. 
         [0005]    Within radiotherapy treatment planning, a search is motivated for a patient-specific best compromise between applying sufficient dose to the target and sparing organs at risk and healthy tissue. Typically, a list of clinical goals (“treatment protocol”) is aspired to, containing (minimally required to optimal) target doses for specific volume quantiles of the planning structures, or reversely target quantiles for specific dose levels. 
         [0006]    In conventional planning, (surrogate) objectives for the different aspects are weighted and combined in a single objective function, and the plan is optimized with respect to this function. However, the best weighting factors for the specific patient are not known a-priori and must be found out by the planner by trial-and-error (“human iteration loop”). 
         [0007]    Therefore it would be desirable to have a system and method that allows treatment planners and physicians to directly and systematically explore the range of treatment options with respect to the criteria given by the treatment protocol. This would allow them to understand the trade-offs for the individual patient while avoiding the time-consuming human iteration loop. 
       SUMMARY OF THE INVENTION 
       [0008]    The invention overcomes the drawbacks of conventional planning (“human iteration loop”). It satisfies the need to calculate well-placed points on a Pareto surface (PS) also in the case that parts or all of the PS are non-convex, as results from the usage of non-convex objective functions or a non-convex feasible set. The invention applies the claimed method to intensity modulated radiation therapy (IMRT) and volumetric arc therapy (VMAT) problems, where it facilitates planning with respect to quantile based criteria. 
         [0009]    The invention includes an algorithm that uses a feasible set and objective functions and this condition is satisfied by staying within a linear programming environment or within the context of convex feasible set and convex objective functions. 
         [0010]    The invention provides for a computer program product. It further provides 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 define adjustments of a radiation apparatus. It further provides a method to select a desired portion of a subject to receive a desired radiation dose. Further the invention provides for several systems that select a desired portion of a subject to receive a (desired) radiation dose. Further a system as according to the invention both, configured and adapted to operate for selecting m objective functions and an initial m-dimensional starting box believed to contain a Pareto surface. Other systems and methods are disclosed as well in the other independent claims. 
         [0011]    The Pareto Surface determines the portions of a subject (a patient) to receive a determined radiation dose. 
         [0012]    The computer program product provides code for being executed on a computer. Preferably a multi-kernel (or multi core) computer is used. The executing code on the computer will satisfy functional tasks and requirements, such as selecting a desired portion of a subject to receive a radiation dose, by determining or developing a previously not known multi-dimensional surface. This surface is shaped by Pareto optimal points that are situated in the objective space. The Pareto points define adjustments of a given radiation apparatus, such it is not only a point but a whole set of physical or radiation parameters that give a specific setting of this radiation apparatus. 
         [0013]    The general approach is to determine a plurality Pareto points on a Pareto surface or the other way around, the plurality of Pareto (optimal) points defines the Pareto surface. Each point however has to be determined, as it is not yet known and when known it is one item of the Pareto surface. 
         [0014]    Therefore, the claim also mentions the assumed Pareto surface which is a surface that is not yet known but only assumed and will become more determined by each additional point that is found and defined to be a proper point of the real Pareto surface. In other terms, the assumed Pareto surface or the un-known Pareto surface will be built according to the invention into a real Pareto surface by adding up many points that span the real Pareto surface. 
         [0015]    Reference points are used therefor and there is a development direction for each reference point that is used. Each (or at least most) specific step of the individual optimization run is commencing with step (b), continued with step (c) and resulting in step (d), where a new determined final front point is found that represents a setting of the radiation apparatus. 
         [0016]    In step (b) a first reference point is selected, as mentioned, on a backside of an assumed Pareto surface (the assumed Pareto front). The first direction is emerging therefrom, towards the assumed Pareto surface from behind. In step (c) a first starting adjustment is provided and an iterative process is used to develop forward a minimum criterion in steps until a final adjustment is reached. Developing forward a minimum criterion is mathematically correct, although starting at p i  and running towards the Pareto surface (the negative of the arrow q i ) from behind. 
         [0017]    The minimum criterion can be the one from Pascoletti-Serafini, see e.g. Scalarizing vector optimization problems by Pascoletti and Serafini, Journal of Optimization Theory and Applications, Vol. 42, 1984, Issue 4, pp. 499 to 524 (herein Pascoletti-Serafini). 
         [0018]    This final adjustment is still able to be implemented on the radiation apparatus, and thus stopping the forward development with the minimum criteria provides for a determination of the achieved final adjustment (developed from the starting adjustment of the radiation apparatus). This final adjustment is represented by the final front point for this individual optimization run (and found to be a real Pareto point of the assumed Pareto surface or Pareto front). 
         [0019]    In step (d) along the direction towards to the front and towards to the back (in front and behind) of the achieved final front point certain portions of the objective space are dismissed. The direction just says that back and forth there is space that cannot contain any other portion that belong to the real Pareto surface. The definition of a Pareto point is that it cannot developed further forward, unless making another parameter adjustment (or setting) less suitable or worse and thus the mathematical definition says the front point must be the Pareto point. Any space towards the back of the direction is also no objective space that may or will contain parts of the Pareto surface. Therefore those spaces can be dismissed or removed or not taken further care of. 
         [0020]    In an exemplary approach, horizontal lines or planes can be used that are building boxes in front and behind of the final front point, each box having a horizontal and vertical plane crossing in the final front point. Thus they may be removed from further consideration as volume of those boxes. 
         [0021]    As the boxes used can be two- or more than two-dimensional, hyperboxing is used as technical term for step by step removing parts from the objective space, to finally retain the Pareto surface out of the previously full objective space. Compared to the volume, this surface is only “thin” along with a stopping criterion that determines the approximation of this surface. The smaller the stopping criterion is the thinner and precise the Pareto front is built. 
         [0022]    Repeating steps (b) to (d) will give more and more Pareto points (PoP) from the previously unknown or previously assumed Pareto surface to build a real Pareto surface. 
         [0023]    Having dismissed certain portions for each individual optimization run, the program or the code continues with remaining, now more determined portions that are assumed to each contain a part of the real Pareto surface presently only being an assumed Pareto surface. 
         [0024]    In an embodiment from previously calculated points on the Pareto surface, the invention finds boxes such that the Pareto surface is contained in their union. The box dimensions provide a measure for the approximation quality realized by the previously calculated points. The box dimensions may also provide the means to determine where a new point on the PS should be calculated. 
         [0025]    These benefits are not available in prior art approximation methods, such as the normal boundary intersection method, the epsilon constraint method, or the uniform weights scalarization approach. 
         [0026]    These benefits are not available for a non-convex PS (non-convex objectives functions or non-convex feasible set) by means of the sandwiching method, see U.S. Pat. No. &#39;316 by Grant et al. mentioned previously. 
         [0027]    In accordance with the present disclosure, a method for selecting a desired portion of a subject to receive a radiation dose is disclosed. The method includes selecting a plurality of points on a Pareto surface (PS) using a plurality of reference points and directions to identify the plurality of points, choosing additional reference points and directions to iteratively run to build up a model of the Pareto surface (or even the surface as Pareto front as such). The method further includes identifying new reference points and directions until a stopping criterion is met. 
         [0028]    In accordance with another aspect of the present invention, a method for selecting a desired portion of a subject to receive a radiation dose is disclosed. The method includes enclosing the Pareto surface in a union of boxes and calculating box dimensions to gauge the approximation quality, and repeating the above steps until a stopping criterion (i.e. criteria) is met. 
         [0029]    In accordance with another aspect of the invention, a system for selecting a desired portion of a subject to receive a radiation dose is disclosed. The system includes a computer having a computer program that when executed causes the computer to select a plurality of points believed to be on a Pareto surface about the desired portion of the subject using a plurality of reference points and directions to identify the plurality of points, iteratively choosing additional reference points and directions to run to build up a model of the Pareto surface. The computer is also caused to identify new reference points and directions (for a multiple kernel several in parallel) and repeat the steps until a stopping criterion is met. 
         [0030]    In another embodiment a computer program product is disclosed, usable on a computer. 
         [0031]    The product comprising a 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 the computer. The executing code is configured for . . . 
         [0000]    a) determining a plurality of Pareto points on a Pareto surface (PS) using a plurality of reference points and development directions, emerging therefrom;
 
b) selecting a first reference point (named p 1  as it is the first) on a back side of an assumed Pareto surface (PS′) and first direction (named q 1  as it is the first) emerging therefrom towards the assumed Pareto surface from behind;
 
c) selecting a first starting adjustment and iteratively developing forward a minimum criterion in steps until a final adjustment is reached that still is implementable in the radiation apparatus; and stopping the forward development, thereby determining the achieved final adjustment represented by a final front point as a real Pareto point of the assumed Pareto surface;
 
d) along the first direction, dismissing undetermined portions of the objective space in front of and behind of the achieved final front point as not containing parts of the Pareto surface, and continuing with other remaining more determined portions that are assumed to each contain a part of the assumed Pareto surface.
 
         [0032]    In another embodiment a method is disclosed. The method determines or develops a previously unknown multi-dimensional surface shaped by Pareto optimal points, situated in the objective space; the Pareto optimal points define adjustments of a radiation apparatus. The method has the following steps. 
       Determining Pareto Points for a Pareto Surface (PS). 
       [0033]    Select a first reference point (named p 1  as it is the first point) on a back side of an assumed, not yet defined Pareto surface (PS′) and selecting a first direction (named q 1  as it is the first direction) emerging therefrom towards the assumed Pareto surface (PS′) from behind. 
         [0034]    Select a first starting adjustment and iteratively developing forward along the first direction a minimum criterion until a final adjustment is reached that still is implementable in the radiation apparatus; and stopping the forward development, thereby determining the achieved final adjustment represented by a final front point as a real Pareto point of the assumed Pareto surface (PS′). 
         [0035]    Dismiss undetermined portions of the objective space in front of and behind of the achieved final front point as not containing parts of the real Pareto surface. 
         [0036]    In another embodiment a System is disclosed. The system is to select a desired portion of a subject to receive a radiation dose. The system comprises . . . 
         [0000]    A computer and a non-transitory computer readable medium having computer readable program code thereon, the computer readable program code comprising a series of computer readable program steps to optimize, with the use of the computer, said radiation dose by steps, at least the following ones . . .
 
a) covering the yet un-determined Pareto surface by a set of m-dimensional boxes, the Pareto surface having a plurality of points wherein at least some points from the plurality of points are selected by solving a Pascoletti-Serafini problem along the diagonal of a box;
 
b) finding the largest box according to some measure;
 
c) comparing the measure to a threshold value;
 
and repeats steps a) through c) until the measure of the largest box is below the threshold value.
 
         [0037]    In still another embodiment a further System is disclosed. The system is for selecting a desired portion of a subject to receive a radiation dose, the system comprising: a computer having a computer program that when executed causes the computer to . . . 
         [0000]    a) select a plurality of points believed to be on a Pareto surface (PS) about the desired portion of the subject using a plurality of reference points and directions, wherein at least some points from the plurality of points are selected by solving a Pascoletti-Serafini problem;
 
b) choose an additional reference point and a corresponding direction to run to build up an additional point of the model of the Pareto surface;
 
c) identify a new reference point and direction including a reference point being the upper corner and the direction being the diagonal of a m-dimensional box containing part of the Pareto surface;
 
and iteratively repeats steps b) and c) using the new reference point and direction for each repetition until a stopping criterion is met.
 
         [0038]    In still another embodiment a further method is disclosed. The method is for selecting a desired portion of a subject to receive a radiation dose. 
         [0039]    A desired portion of a subject is selected, to receive a radiation dose. A computer system is provided, and with the use of the computer system, data, associated with the radiation dose, is received and the received data is processed by . . . 
         [0000]    a) selecting surface points of a Pareto surface (PS) that defines m objective functions, and an initial m-dimensional starting box containing the points;
 
b) computing, from the initial box and points, a set of outer knee points and inner knee points determining a set of m-dimensional boxes whose union covers the Pareto surface;
 
c) increasing a number of points of the Pareto surface by adding a point in one of the boxes whose union covers the Pareto surface.
 
     
    
     
       INTRODUCTION TO THE DRAWINGS 
         [0040]    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. 
           [0041]      FIG. 1  shows an example of a workflow of the disclosed method. 
           [0042]      FIG. 1 a    shows a radiation apparatus  60 ,  62  and  69 , taken from U.S. Pat. No. 7,391,026. 
           [0043]      FIG. 2  exemplifies how radiation dose affects multiple volumes of interest (VOIs). With each VOI a dose evaluating objective function ƒ i  is associated, measuring some dose average. 
           [0044]      FIGS. 3 a  to 3 d    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  51   a ,  51   b  in  FIG. 3 b   ;  5 ,  5  and  52   a ,  52   b  in  FIG. 3   d.    
           [0045]      FIGS. 3 e  to 3 h    further explain  FIG. 3 d    with respect to inner knee points and outer knee points as well as their sets and the sets of boxes B 1 , B 2 , B 3  and B 4  that knee points define. The initial box used for starting is shown stippled in  FIG. 3 h    (and  FIGS. 3 f  and 3 g   ). 
           [0046]      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. 
           [0047]      FIG. 5  display Pareto surface approximations, obtained by hyperboxing . . . 
           [0048]      FIG. 5 a    shows the points with the ideal point (0,0,0) in front 
           [0049]      FIG. 5 b    shows the triangulated points seen from a side. 
           [0050]      FIG. 6  are alternative Pareto surface examples of triangulated Pareto surfaces PS as developed by examples of the disclosed method.  FIGS. 6 a  to 6 e    show the PS from different sides. 
           [0051]      FIG. 7  shows a computer with multiple cores or kernels used in the methods disclosed herein. 
       
    
    
     DETAILED DESCRIPTION OF THE EMBODIMENTS 
       [0052]    One approach is the planning problem in intensity modulated radiation therapy (IMRT) and in volumetric arc therapy (VMAT). 
         [0053]    The radiation apparatus is not shown, as suitable ones can be found in the prior art for IMRT and VMAT therapy. 
         [0054]    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 
         [0000]      min xεX {ƒ i ( x )}, i= 1, . . . , m  
 
         [0000]    where
 
x is the decision variable
 
X denotes the feasible set
 
ƒ i (x) are the objectives such that smaller values denote preferable properties of x.
 
         [0055]    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”. 
         [0056]    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. 
         [0057]    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). 
         [0058]    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). 
         [0059]    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&#39;s preferences beforehand, such as solving weighted sum problems with fixed weights, goal programming, and lexicographic ordering. 
         [0060]    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. 
         [0061]    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. 
         [0062]    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. 
         [0063]    From the representation points, a surface approximation of the Pareto surface PS can be obtained by triangulation schemes, such as the Delaunay triangulation. 
         [0064]    Individual solutions to the multi-objective problem can be computed from scalarizations. 
         [0065]    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. 
         [0066]    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 
         [0000]      α→min
 
         [0000]        s.t. p+αq=F ( x )+λ
 
         [0000]      (α, x ,λ)ε         × X×             +   m 
 
         [0000]      Where . . . 
         [0000]        F ( x )=( f   1 ( x ), f   2 ( x ), . . .  f   m ( x )) 
         [0000]    is the vector of objective function values and p and q are the parameters. 
         [0067]    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. 
         [0068]    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&#39;s 95% quantile. Translating these clinical goals into objectives yields non-convex functions and a non-convex PS. 
         [0069]    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. 
         [0070]    Because of this usage of m-dimensional boxes, we refer to the algorithm as “hyperboxing”. 
         [0071]    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. 3   a.    
         [0072]    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. 
         [0073]    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. 
         [0074]    Applied to  FIG. 3 b   , the following iteration chooses the largest box  51   a  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 1  of the chosen box is solved, the solution is added as a point y 11  to the representation set, and the set of covering boxes is updated again. This applies in the next iteration with y 31  and y 21  in  FIG. 3 c   , where the diagonals are q 2  and q 3 . The iteration using the cut away boxes  5  and  52   a  is shown in  FIG. 3 d   . Points y 31  and y 21  are added to the representation set. The set of covering boxes is then updated again. 
         [0075]    This iteration is repeated, until a stopping criterion (or stop criterion) with respect to a size of the largest box is met. 
         [0076]    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. 3 a    to  FIG. 3 d   , as mentioned before. 
         [0077]    Assume having obtained (α, x 11 , λ) as the solution of the Pascoletti-Serafini problem along the diagonal q 1  of an initial box, as depicted in  FIG. 3 a   . This defines the Pareto point y 11 =F (x 11 ). It holds that there are no Pareto optimal points in the area left-hand and below the point y 11  (the point marked square, not round on the diagonal q in  FIG. 3 b   ). Otherwise a lower value for α 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). 
         [0078]    Also, there are no Pareto optimal points in right-hand and above y 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. 3   b.    
         [0079]    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. 3 c    and  FIG. 3   d.    
         [0080]    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. 
         [0081]    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. 
         [0082]    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”. 
         [0083]    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 . 
         [0084]    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. 
         [0085]      FIG. 1  shows one preferred workflow of hyperboxing.  FIG. 1 a    shows a radiation apparatus as taken e.g. from prior art U.S. Pat. No. 7,391,026 (Fraunhofer) and its  FIG. 1 . 
         [0086]    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. 
         [0087]    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. 
         [0088]    While it is possible that the Pareto surface PS exceeds this region, it is a reasonable estimate to start with. 
         [0089]    Finding the true minimal starting box that contains the PS is a computational problem in dimensions greater than two. 
         [0090]    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 . . . 
         [0000]      μ( B )=min{ L   i   ,i= 1, . . . , m}.  
 
         [0091]    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. 
         [0092]    To find the point, the Pascoletti-Serafini-Problem along q i  of the chosen box is solved. “i” is the iteration index of loop of steps  102  to  105 . 
         [0093]    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. 3   c.    
         [0094]    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.). 
         [0095]    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. 
         [0096]    The method then proceeds to  102  again with determining the largest box from the updated set of boxes at  105 . 
         [0097]    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. 
         [0098]    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. 
         [0099]    This allows parallelization on a computer  80  with multiple cores (or multiple kernels) as shown in  FIG. 7  and  FIG. 1 a   . 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 1 , q 2 , q 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 . 
         [0100]    Applying  FIG. 4  to  FIG. 7  (or vice versa), calculation along q 1  could be performed in core  1 , calculation along q 2  could be performed in core  2  and calculation for direction q 3  could be performed in core  3 . Thus at least a few calculations along a few directions q 1 , q 2  and q 3  are associated with separate kernels of the shown multi-kernel processor  80  to provide parallel processing. 
         [0101]    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. 
         [0102]    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. 
         [0103]    The computerized system  80  is schematically shown in  FIG. 1 a    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. 1 a    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. 
         [0104]    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. 
         [0105]    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. 
         [0106]    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 U.S. Pat. No. &#39;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. 
         [0107]    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 . 
         [0108]    When selecting a desired portion of a subject to receive a radiation dose provided by the beam head  64  in  FIG. 1 a   ,  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 f i , or is associated with it, measuring a dose average for the respective volume of interest. Thus the objective function f 1  is related to a dose average f 1 . The same applies for volume of interest  2  which is associated with the objective function f 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 f 3  and f 4 . 
         [0109]    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. 
         [0110]    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. 
         [0111]    It may be specified, how to add a point to the Pareto surface PS to best improve the representation of the PS. 
         [0112]    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. 
         [0113]    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. 
         [0114]    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. 3 c   . 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. 3 c    and  FIG. 3   d.    
         [0115]    This may, for example, be achieved by updating both, the set of so-called “outer knee points” v i  and the set of so-called “inner knee points” w i  as present in  FIG. 3 d    and highlighted in  FIGS. 3 f  and 3 g   , both together in  FIGS. 3 e    and  3   h.    
         [0116]    The outer knee points v of  FIG. 3 e    form a set of points v 1 , v 2 , v 3 , v 4 , . . . with the following property: for each point (e. g. y 11 , y 21 , y 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. 3 e    form a set with the following property: for each point (e. g. y 11 , y 21 , y 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. 
         [0117]    The updated set of boxes is then defined by the set of pairs comprising one inner knee point w i  and one outer knee point v j  from the updated sets, respectively, such that the two points span a box with non-zero volume. 
         [0118]    When viewing at  FIG. 3 e   , the outer knee points v and the inner knee points w of  FIG. 3 d    can be easily spotted. In  FIG. 3 e    they have names, v j  (j is running from 1 to 4) and w i  (i is running from 1 to 4). In  FIG. 3 d    the names of the knee points are not yet presented. In both  FIG. 3 e    and  FIG. 3 h    the boxes are each spanned by a pair of knee points. Namely, these are box B 1 , spanned by knee points v 1  and w 2  (one outer and one inner knee point), box B 2  is spanned by knee points v 2  and w 2 , box B 3  is spanned by knee points v 3  and w 3  and box B 4  is spanned by knee points v 4  and w 4 . These are boxes of non-zero volume that in their entirety form a set of covering boxes that approximate the Pareto surface PS. 
         [0119]    Now separating  FIG. 3 e    (or  FIG. 3 d   ) will yield  FIGS. 3 f  and 3 g   .  FIG. 3 g    shows the inner knee points w 1 , w 2 , w 3 , w 4  and their relation to a respective point on the Pareto surface.  FIG. 3 f    shows the outer knee points v 1 , v 2 , v 3 , v 4 , also related to a respective point on the Pareto surface. The points on the Pareto surface are 
         [0120]    Using  FIGS. 3 g  and 3 f    together and putting them on top of each other, will yield  FIG. 3 h   , where the same structure is given as shown in  FIG. 3 d   , outlining the boxes covering the Pareto front PS, as discussed before. This is exactly what is shown in  FIG. 3 d   , just outlining or enhancing the points that are present in  FIG. 3 d    already. 
         [0121]    From  FIGS. 3 d , 3 e    it may be taken that from the initial box B 0  (given by broken line in  FIG. 3 h   ) and the set of Pareto points y 11 , y 21  and y 31  a set of outer knee points v 1 , v 2 , v 3 , v 4  and a set of inner knee points w 1 , w 2 , w 3 , w 4  may be computed. 
         [0122]    The initial box B 0  is shown in broken lines in  FIG. 3 h    (encompassing the set of boxes B 1  to B 4 , 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 B 1  to B 4  boxes whose union covers the Pareto surface. 
         [0123]    Thus surface points on the Pareto surface PS define m objective functions, using the initial m-dimensional starting box containing the points. 
         [0124]    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. 
         [0125]    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&#39;s reach (PHOSITO&#39;s reach) and thus within the scope of the invention as claimed. Therefore, the disclosure should not be limited to a particular embodiment.