Abstract:
A method of multileaf collimator (MLC) leaf positioning in tracking-based adaptive radiotherapy is provided. The method includes determining a radiotherapy beam pattern by transforming a treatment beam plan into radiotherapy beam coordinates, determining a dose discrepancy between the radiotherapy beam pattern and a deliverable MLC aperture, where the dose discrepancy includes a sum of an overdose cost and an underdose cost to a treatment volume, and minimizing the dose discrepancy, where the dose discrepancy minimization provides a determined deliverable MLC aperture for the radiotherapy beam.

Description:
FIELD OF THE INVENTION 
       [0001]    The current invention relates to radiotherapy. More particularly, the invention relates to systematically adjust a multileaf collimator (MLC) to conform to the geometric and topological motion of a treatment volume encompassing the target during the course of radiotherapy. 
       BACKGROUND OF THE INVENTION 
       [0002]    Radiotherapy aims to deliver an ablative dose to the tumour with minimal normal tissue exposure. However, anatomical motion during treatment results in misalignment between beam and target, compromising treatment efficacy. An ideal motion management strategy is to reposition and reshape the beam aperture in response to instantaneous motion. There exists a need to provide an optimal deliverable beam pattern. For multileaf collimator (MLC) tracking based delivery, a systematic leaf adaptation method is desired. 
         [0003]    Given an estimate of instantaneous motion, it is desirable to transform the planning beam according to the such motion to generate an ideal aperture. However, such an ideal aperture is often undeliverable, due to physical constraints such as finite MLC leaf widths, and more importantly, the paired leaf structure. 
         [0004]    What is needed is an optimization framework that provides a deliverable MLC configuration that is closest to the ideal aperture, where closeness is defined rigorously as the cumulative cost in terms of underdose to target and overdose to healthy tissue. 
       SUMMARY OF THE INVENTION 
       [0005]    To address the shortcomings in the art, a method of multileaf collimator (MLC) leaf positioning in tracking-based adaptive radiotherapy is provided. The method includes determining a radiotherapy beam pattern by transforming a treatment beam plan into radiotherapy beam coordinates, determining a dose discrepancy between the radiotherapy beam pattern and a deliverable MLC aperture, where the dose discrepancy includes a sum of an overdose cost and an underdose cost to a treatment volume, and minimizing the dose discrepancy, where the dose discrepancy minimization provides a determined deliverable MLC aperture for the radiotherapy beam. 
         [0006]    In one aspect of the invention, the radiotherapy beam coordinates are based on the projection of the translation, rotation and/or deformation of all or part of a radiotherapy beam pattern into the beam coordinates. 
         [0007]    In another aspect of the invention, the overdose cost includes an integration of pixelwise overdose penalties, where the overdose penalties can be spatially variant. 
         [0008]    In a further aspect of the invention, the underdose cost includes an integration of pixelwise underdose penalties, where the underdose penalties can be spatially variant. 
         [0009]    In yet another aspect of the invention, at least one previously determined dose discrepancy is used when determining a next deliverable MLC aperture for the radiotherapy beam. 
         [0010]    In a further aspect of the invention, a determination of the overdose cost and a determination of the underdose cost are based on a tissue type and a radio-sensitivity of the tissue. 
     
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         [0011]      FIGS. 1   a - 1   b  show an undeliverable beam aperture. 
           [0012]      FIG. 2  shows a schematic of a leaf fitting configuration, according to the current invention. 
           [0013]      FIGS. 3   a - 3   d  show automatic leaf sequencing with the optimization method, according to the current invention. 
           [0014]      FIGS. 4   a - 4   h  show an original and undeliverable IMRT field and the deliverable beam under various overdose/underdose tradeoff configurations, according to the current invention. 
       
    
    
     DETAILED DESCRIPTION 
       [0015]    Radiation therapy aims to deliver an ablative radiation dose to tumor targets while sparing the surrounding normal healthy tissue. While current development in radiotherapy machines enables focused radiation beams to be directed with high precision, anatomical motion during treatment causes misalignment between the external beam and the internal anatomy, compromising the treatment efficacy. Therefore, it is important to reduce the impact of motion in radiotherapy. To achieve this goal, one aspect of the current invention synchronizes the treatment with anatomical motion, broadly known as “tracking”. 
         [0016]    The key idea behind tracking-based adaptive radiotherapy is to modify the treatment plan and reposition the beam to compensate for the anatomical motion. When rigid translational anatomical motion is assumed, there is a simple scheme to move the MLC to “follow” the target. When general motion is considered, however, a more sophisticated and systematic scheme is required for MLC adaptation. The most ideal situation would be performing a full-fledged re-planning in real-time. However, this is impractical due to heavy computation and the as-yet unresolved quality assurance issue. The problem of leaf arrangement for real-time intensity modulation has been studied quite intensively. The current invention applies a first order approximation to this ideal situation by adapting the MLC leaf configuration derived from a treatment plan to the instantaneous motion. The planned beam aperture is transformed with a collapsed transformation describing online motion on the beam-eye-view (BEV) plane to generate a new aperture. However, motion could transform a deliverable beam aperture into quite complex shapes that are no longer deliverable, due to paired leaf structure with finite leaf resolution. More specifically, this occurs when a single leaf track intersects with multiple open apertures. Referring now to the figures,  FIGS. 1   a - 1   b  show a schematic drawing of an MLC system  100 , where multiple leafs  102  are configured to provide at least one aperture  104  that is an optimization-based leaf sequencing methodology that outputs a deliverable leaf configuration that best approximates the ideal pattern.  FIG. 1   a  shows an ideal aperture  104  is generated by translating a deliverable planned aperture by a fraction of the width of a leaf  102 . The single leaf track  106  in  FIG. 1   a  is enlarged in  FIG. 1   b , which intersects two open apertures  104 . There exists no placement of the leaf  102  that could exactly realize this pattern. In fact, with the motion perpendicular to direction of the leaf  102  being a fraction of the width of the leaf  102 , no exact placement of the leaf  102  exists for any of the leaves  102  intersecting ideal aperture  104 . 
         [0017]    Another challenge for tracking-based adaptive radiotherapy is the differential motion between tumor target and critical structures. The majority of existing studies for tracking based treatment adaptation are based on the assumption that the complete region of interest (ROI) undergoes similar motion, so the relative anatomical configuration is maintained throughout treatment. The feasibility of this assumption very much depends on the site, where it is more valid for prostate but is questionable for treating thoracic tumors, which may move towards the spinal cord with respiration. Therefore, it is desirable to have a leaf adaptation method that automatically accounts for the potential varying structural characteristics of the ROI. This property is embedded naturally in the optimization aspect of the current invention. 
         [0018]    The current invention addresses the leaf sequencing problem described above in the presence of motion. According to one aspect, a plan aperture f and an estimated anatomical motion T are collapsed onto the BEV following the methodology described above. The ideal motion-compensated aperture is given by composing the plan aperture with the collapsed BEV motion g=f·T. According to the invention, this map is represented with a binary function over the ROI: Ω→{0,1}, so that 
         [0000]    
       
         
           
             
               g 
                
               
                 ( 
                 
                   x 
                   _ 
                 
                 ) 
               
             
             = 
             
               { 
               
                 
                   
                     1 
                   
                   
                     
                       
                         
                           x 
                           _ 
                         
                         ∈ 
                         
                           transformed 
                            
                           
                               
                           
                            
                           plan 
                            
                           
                               
                           
                            
                           aperture 
                            
                           
                               
                           
                            
                           opening 
                         
                       
                       ; 
                     
                   
                 
                 
                   
                     0 
                   
                   
                     
                       else 
                       , 
                     
                   
                 
               
             
           
         
       
     
         [0000]    where  x =(x, y) denotes the beam element location in the BEV. Since leaf sequencing is of interest to approximate the desired beam pattern g, it suffices to consider the ROI in the 2-dimensional BEV coordinate. 
         [0019]    Without loss of generality, x-coordinate is aligned with the leaf tracking direction. Let Δ be the leaf resolution along y-direction and the complete aperture is covered by the range of N leaves, starting from y=0. The problem is restated as: finding the best N pairs of leading and trailing locations that delivers an aperture that is closest to the ideal beam pattern given by g. The leaf positions are parameterized as {α i   L , α i   T } where the subindex i=1, 2, . . . , N indexes the leaf pair and the superscripts L and T indicate the leading and trailing leaf respectively. 
         [0020]    To rigorously define an optimization objective, the “dose” discrepancy between a deliverable aperture corresponding to {α i   L , α i   T } i=1   N  the desired beam pattern g are qualified. Such a discrepancy is characterized as the sum of underdose and overdose cost. A pair of unitary underdose and overdose costs λu( x ), λo( x ) are associated with each voxel location, based on its tissue type and radio-sensitivity. If  x  belongs to tumor region, then the local underdose cost λu( x ) should be high and the local overdose cost λo( x ) should be low. Conversely, for radiosensitive healthy tissues that should be prevented from radiation, λu( x ) should be low and λo( x ) high. A tumor voxel desired to be irradiated (g( x )=1) but blocked by the delivered MLC pattern contributes an underdose cost λu( x ) to the overall discrepancy; conversely, a healthy tissue voxel desired to be shielded (g( x )=0) but exposed by the delivered MLC pattern contributes an overdose cost λo( x ) to the overall discrepancy.  FIG. 2  shows a schematic of the leaf fitting setup  200 . 
         [0021]    Assume the leaf tracks are of widths Δ, then the underdose cost for leaf track i is given by: 
         [0000]      Φ i,u =∫   x :xε[L,α     i       L     )∪(α     i       T     ,R],yε[(i−1)Δ,iΔ λ u (   x   ) g (   x   );
 
         [0022]    and the overdose cost reads: 
         [0000]      Φ i,o =∫   x :xε(α     i       T     ,α     i       T     ),yε[(i−1)Δ,iΔ λ o (   x   )(1− g (   x   );
 
         [0023]    Where L and R indicate the left and right jaw location. 
         [0024]    To reveal the dependence of Φi,o and Φi,u on the parameters {α i   L ,α i   T }, the indicator function I is invoked: 
         [0000]    
       
         
           
             
               I 
                
               
                 ( 
                 u 
                 ) 
               
             
             = 
             
               { 
               
                 
                   
                     1 
                   
                   
                     
                       
                         u 
                         &gt; 
                         0 
                       
                       ; 
                     
                   
                 
                 
                   
                     0 
                   
                   
                     
                       else 
                       . 
                     
                   
                 
               
             
           
         
       
     
         [0025]    Then the underdose and overdose cost can be rewritten as 
         [0000]      Φ i,u =∫   x :yε[(i−1)Δ,iΔ) λ u (   x   ) g (   x   )[ I (α i   L   −x )+ I ( x−α   i   T )];
 
         [0000]      Φ i,o =∫   x :yε[(i−1)Δ,iΔ) λ o (   x   )(1− g (   x   ) I ( x−α   i   L )(α i   T   −x ).
 
         [0026]    The objective is to find {α i   L , α i   T } to minimize the overall objective: 
         [0000]    
       
         
           
             Φ 
             = 
             
               
                 
                   ∑ 
                   i 
                   
                       
                   
                 
                  
                 
                     
                 
                  
                 
                   Φ 
                   
                     i 
                     , 
                     u 
                   
                 
               
               + 
               
                 
                   Φ 
                   
                     i 
                     , 
                     o 
                   
                 
                 . 
               
             
           
         
       
     
         [0027]    For simplicity, the leaf velocity condition is ignored in the introduction of basic formulation, where this condition can be incorporated by imposing constraints on the optimization problem. 
         [0028]    In the absence of any additional constraint, the summation form of the objective function enables decoupled optimization for each leaf pair, i.e., it suffices to minimize 
         [0000]      Φ i =Φ i,u +Φ i,o .
 
         [0029]    From here, a specific leaf track is studied and all integrations are restricted to this range. Note that 
         [0000]      Φ i,u =∫   x :yε[(i−1)Δ,iΔ) λ u (   x   ) g (   x   )[ I (α i   L   −x )+ I ( x−α   i   T )];
 
         [0000]      Φ i,o =∫   x :yε[(i−1)Δ,iΔ) λ o (   x   )(1− g (   x   ) I ( x−α   i   L )(α i   T   −x ).
 
         [0030]    indicate that the integrand depends on the y coordinate via λ and g, which is independent of the parameters {α i   L , α i   T }. This observation allows one to first integrate along the y-direction, and reduces the problem to a 1-dimensional optimization. More specifically, for each leaf track i, one can define the corresponding 1D function 
         [0000]        c   i,u ( x )=∫ (i−1)Δ   iΔ λ u (   x   ) g (   x   ) dy  
 
         [0000]        c   i,o ( x )=∫ (i−1)Δ   iΔ λ o (   x   )(1− g (   x   ) dy.  
 
         [0031]    These definitions reduce the cost function to 
         [0000]      Φ u,i   =∫c   i,u ( x )[( I (α i   L   −x )+ I ( x−α   i   T )] dx  
 
         [0000]      Φ o,i   =∫c   i,o ( x ) I ( x−α   i   L ) I (α i   T   −x ).
 
         [0032]    It can be shown with simple arithmetic manipulation that minimizing Φ i  is equivalent to minimizing 
         [0000]      Ψ i =∫ α     i       L     α     i       T     [−c   i,u ( x )+ c   i,o ( x )] dx.  
 
         [0033]    Let 
         [0000]    
       
         
           
             
               
                 c 
                 i 
               
                
               
                 = 
                 Δ 
               
                
               
                 ( 
                 
                   
                     c 
                     
                       i 
                       , 
                       u 
                     
                   
                   + 
                   
                     c 
                     
                       i 
                       , 
                       o 
                     
                   
                 
                 ) 
               
             
             , 
           
         
       
     
         [0000]    then the first order necessary condition for optimality is given by 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       ∂ 
                       
                           
                       
                     
                     
                       ∂ 
                       
                         α 
                         i 
                         L 
                       
                     
                   
                    
                   
                     Ψ 
                     i 
                   
                 
                  
                 
                    
                   
                     
                       α 
                       i 
                       L 
                     
                     = 
                     x 
                   
                 
               
               = 
               
                 
                   
                     c 
                     i 
                   
                    
                   
                     ( 
                     x 
                     ) 
                   
                 
                 = 
                 0 
               
             
             ; 
           
         
       
       
         
           
             
               
                 
                   
                     ∂ 
                     
                         
                     
                   
                   
                     ∂ 
                     
                       α 
                       i 
                       T 
                     
                   
                 
                  
                 
                   Ψ 
                   i 
                 
               
                
               
                  
                 
                   
                     α 
                     i 
                     T 
                   
                   = 
                   x 
                 
               
             
             = 
             
               
                 
                   - 
                   
                     
                       c 
                       i 
                     
                      
                     
                       ( 
                       x 
                       ) 
                     
                   
                 
                  
                 
                   ( 
                   x 
                   ) 
                 
               
               = 
               0. 
             
           
         
       
     
         [0034]    In addition, Φ i  needs to be locally convex at the optimal {α i   L , α i   T }. Fortunately, the leading and trailing leaf positions are only related by the constraint that (α i   L &lt;α i   T ), and otherwise decoupled. This implies a diagonal Hessian matrix when Φ i  is twice differentiable, and the positive-definiteness of the Hessian reduces to the positiveness of the diagonal elements 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       ∂ 
                       2 
                     
                      
                     
                         
                     
                   
                   
                     ∂ 
                     
                       
                         ( 
                         
                           α 
                           i 
                           L 
                         
                         ) 
                       
                       2 
                     
                   
                 
                  
                 
                   Ψ 
                   i 
                 
               
               &gt; 
               0 
             
             ; 
           
         
       
       
         
           
             
               
                 
                   
                     ∂ 
                     2 
                   
                    
                   
                       
                   
                 
                 
                   ∂ 
                   
                     
                       ( 
                       
                         α 
                         i 
                         T 
                       
                       ) 
                     
                     2 
                   
                 
               
                
               
                 Ψ 
                 i 
               
             
             &gt; 
             0. 
           
         
       
     
         [0035]    Substituting the expressions from 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       ∂ 
                       
                           
                       
                     
                     
                       ∂ 
                       
                         α 
                         i 
                         L 
                       
                     
                   
                    
                   
                     Ψ 
                     i 
                   
                 
                  
                 
                    
                   
                     
                       α 
                       i 
                       L 
                     
                     = 
                     x 
                   
                 
               
               = 
               
                 
                   
                     c 
                     i 
                   
                    
                   
                     ( 
                     x 
                     ) 
                   
                 
                 = 
                 0 
               
             
             ; 
           
         
       
       
         
           
             
               
                 
                   
                     
                       ∂ 
                       
                           
                       
                     
                     
                       ∂ 
                       
                         α 
                         i 
                         T 
                       
                     
                   
                    
                   
                     Ψ 
                     i 
                   
                 
                  
                 
                    
                   
                     
                       α 
                       i 
                       T 
                     
                     = 
                     x 
                   
                 
               
               = 
               
                 
                   
                     - 
                     
                       
                         c 
                         i 
                       
                        
                       
                         ( 
                         x 
                         ) 
                       
                     
                   
                    
                   
                     ( 
                     x 
                     ) 
                   
                 
                 = 
                 0. 
               
             
             , 
           
         
       
     
         [0000]    the second order condition reads 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       
                         ∂ 
                         2 
                       
                        
                       
                           
                       
                     
                     
                       ∂ 
                       
                         
                           ( 
                           
                             α 
                             i 
                             L 
                           
                           ) 
                         
                         2 
                       
                     
                   
                    
                   
                     Ψ 
                     i 
                   
                 
                  
                 
                    
                   
                     
                       α 
                       i 
                       L 
                     
                     = 
                     x 
                   
                 
               
               = 
               
                 
                   
                     c 
                     i 
                     ′ 
                   
                    
                   
                     ( 
                     x 
                     ) 
                   
                 
                 &gt; 
                 0 
               
             
             ; 
           
         
       
       
         
           
             
               
                 
                   
                     
                       ∂ 
                       2 
                     
                      
                     
                         
                     
                   
                   
                     ∂ 
                     
                       
                         ( 
                         
                           α 
                           i 
                           T 
                         
                         ) 
                       
                       2 
                     
                   
                 
                  
                 
                   Ψ 
                   i 
                 
               
                
               
                  
                 
                   
                     α 
                     i 
                     T 
                   
                   = 
                   x 
                 
               
             
             = 
             
               
                 - 
                 
                   
                     c 
                     i 
                     ′ 
                   
                    
                   
                     ( 
                     x 
                     ) 
                   
                 
               
               &gt; 
               0. 
             
           
         
       
     
         [0036]    These analyses induce some observations for the behavior of the 1D function c i  which are interpreted as the competing strength of local underdose and overdose cost. These observations include:
       Both α i   L  and α i   T  fall on the zero-crossings of c i .   c i  moves upwards at α i   L , indicating that underdose cost starts to dominate overdose at the leading leaf position. This is expected as the transition from the normal tissue region to the target region occurs.   c i  moves downwards at α i   T , indicating that overdose cost starts to dominate underdose at the trailing leaf position. This marks the transition from the target region to the normal tissue region along the leaf track.   A general observation from combining the above is that tumors with large local underdose cost drive the aperture open to ensure its coverage, and normal tissue regions with high overdose cost shall be shielded with the leaves.       
 
         [0041]    A similar perturbation principle applies to cases when c is nondifferentiable, and yields the following results: 
         [0000]        c (α i   L −)&lt;0 , c (α i   L +)&gt;0;
 
         [0000]        c (α i   L −)&gt;0 , c (α i   T +)&lt;0,
 
         [0000]    where − and + indicate a small perturbation in the negative and positive directions respectively. The zero-crossing conditions are implied by the continuity of the function c and thus dropped for simplicity. A semi-exhaustive algorithm for general inhomogeneous cost is presented in Algorithm 1 below: 
         [0000]    
       
         
               
             
               
               
             
               
               
             
               
               
             
               
               
             
               
               
             
               
               
             
               
               
             
               
               
             
               
               
             
           
               
                   
               
               
                 Algorithm 1 Semi-exhaustive leaf sequencing 
               
               
                 algorithm for general inhomogeneous cost 
               
               
                   
               
             
             
               
                   
               
             
          
           
               
                   
                 for i = 1 to N do 
               
             
          
           
               
                   
                 Compute ID functions c i,u , c i,o  (4) and c i  = c i,u  − c i,o . 
               
               
                   
                 Find set C i   L  = {x : I(c(x + δ))I(−c(x − δ)) = 1}. 
               
               
                   
                 Find set C i   T  = {x : I(c(x + δ))I(−c(x − δ)) = 1}. 
               
               
                   
                 Initialize {circumflex over (α)} i   L  = 0, {circumflex over (α)} i   T  = 0; {circumflex over (Φ)} i  = Φ i (α i   L  , α i   T ). 
               
               
                   
                 for α i   L  ∈ C i   L  do 
               
             
          
           
               
                   
                 for α i   T  ∈ C i   T  and α i   T  &gt; α i   L  do 
               
             
          
           
               
                   
                 {tilde over (Φ)} i  = Φ i (α i   L  , α i   T ). 
               
               
                   
                 if {tilde over (Φ)} i  &lt; {circumflex over (Φ)} i  then 
               
             
          
           
               
                   
                 {circumflex over (α)} i   L  = α i   L , {circumflex over (α)} i   T  = α i   T  ; {circumflex over (Φ)} i  = {tilde over (Φ)} i   
               
             
          
           
               
                   
                 end if 
               
             
          
           
               
                   
                 end for 
               
             
          
           
               
                   
                 end for 
               
             
          
           
               
                   
                 end for 
               
               
                   
                 {{circumflex over (α)} i   L  , {circumflex over (α)} i   T } i=1   N  provides optimal leaf positions. 
               
               
                   
                   
               
             
          
         
       
     
         [0042]    As a special case, homogeneous unitary underdose and overdose cost may be assigned across the ROI. This occurs when insufficient structural information can be used for a localized assignment. The homogeneous assumption, λu( x )=λu and λo( x )=λo, can be used to further reduce the complexity of the general algorithm by pulling the constants out of the integrations, resulting in Algorithm 2 below: 
         [0000]    
       
         
               
             
               
               
             
               
               
             
               
               
             
               
               
             
               
               
             
               
               
             
               
               
             
               
               
             
               
               
             
           
               
                   
               
               
                 Algorithm 2 Semi-exhaustive leaf sequencing 
               
               
                 algorithm for homogeneous unitary cost 
               
               
                   
               
             
             
               
                   
               
             
          
           
               
                   
                 Compute κ = λ u /(λ o  + λ u ). 
               
               
                   
                 for i = 1 to N do 
               
             
          
           
               
                   
                 Compute 1D functions l i (x) = ∫ (i−1)Δ   iΔ  g( x )dy. 
               
               
                   
                 C i  = {x : l i (x) = κΔ}. 
               
               
                   
                 Initialize {circumflex over (α)} i   L  = 0, {circumflex over (α)} i   T  = 0; {circumflex over (Φ)} i  = Φ i (α i   L  , α i   T ). 
               
               
                   
                 for α i   L  ∈ C i  do 
               
             
          
           
               
                   
                 for α i   T  ∈ C i  and α i   T  &gt; α i   L  do 
               
             
          
           
               
                   
                 {tilde over (Φ)} i  = Φ i (α i   L  , α i   T ). 
               
               
                   
                 if {tilde over (Φ)} i  &lt; {circumflex over (Φ)} i  then 
               
             
          
           
               
                   
                 {circumflex over (α)} i   L  = α i   L , {circumflex over (α)} i   T  = α i   T  ; {circumflex over (Φ)} i  = {tilde over (Φ)} i   
               
             
          
           
               
                   
                 end if 
               
             
          
           
               
                   
                 end for 
               
             
          
           
               
                   
                 end for 
               
             
          
           
               
                   
                 end for 
               
               
                   
                 {{circumflex over (α)} i   L  , {circumflex over (α)} i   T } i=1   N  provides optimal leaf positions. 
               
               
                   
                   
               
             
          
         
       
     
         [0043]    To demonstrate the behavior of the optimization method according to the current invention, two illustrative examples are presented, one that demonstrates the effect of inhomogeneous underdose and overdose cost on leaf sequencing; and another presents the adaptation of an IMRT plan when the ROI undergoes rotational motion. 
         [0044]    In this first example, the behavior of the leaf sequencing method is demonstrated in the presence of differential motion between tumor target and surrounding critical structures. An anatomy having a spherical tumor target and an arc-shaped critical structure is simulated. A reference plan is generated when the tumor target is relatively far from the critical structure. Inhomogeneous underdose and overdose unitary cost is assigned according to the discussion above, with a high overdose cost assigned to the critical structure. 
         [0045]    During the simulated treatment process, when the tumour target approaches the critical structure, both the reference plan and the monitored tumour target motion are streamlined into the leaf sequencing algorithm. At each time instant, the ideal aperture is generated as the composition of the initial plan with the target motion.  FIGS. 3   a - 3   d  show snapshots of the evolution of a simulated anatomy  300  having a tumor  302 , and automatically generated leaf patterns  304  for a case with a high critical structure  306  overdose cost. This example mimics the differential motion between lung tumors and the adjacent spinal cord. 
         [0046]    Rotations exceeding 45 degrees for lung tumors and 30 degrees for prostate tumors severely compromise dose delivery. When IMRT treatment is adapted to ROI rotations of such magnitude, the composite ideal aperture is highly likely to be nonconvex along certain leaf tracks, making the beam pattern undeliverable and necessitates a tradeoff between tumor underdose and normal tissue exposure. To illustrate this scenario, assume the only input to the tracking system is the planned IMRT aperture and global rotational information, which can be obtained from electromagnetic tracking or real-time imaging.  FIGS. 4   a - 4   b  shows the original IMRT plan aperture  400  and an ideal aperture shape  402  corresponding to 90 degree rotation, where the choice of 90 degree rotation angle is merely for the convenience of visualization. Since the rotated aperture  402  shape cannot be delivered by the MLC, various underdose/overdose tradeoff results in different aperture patterns that are deliverable. The parameter κ=λo/(λo+λu) describes such a tradeoff. As the unitary underdose cost increases relative to the overdose cost, the optimal aperture pattern tends to be more tolerant towards overdosing the normal regions to ensure target coverage. As shown,  FIGS. 4   c - 4   e  are the optimal deliverable apertures with corresponding κ parameter, and  FIGS. 4   f - 4   h  are the optimal deliverable aperture overlays with the ideal rotated apertures. Further shown are the corresponding overdose  404  and underdose  406  regions. 
         [0047]    The optimization principles of the current invention provide a systematic method to understand and make explicit decisions about the tradeoffs between tumor underdose and normal tissue overdose, when the ideal beam pattern is undeliverable. When the ideal pattern is deliverable, the solution to the optimization problem coincides with it regardless of the specific assignment of tradeoff. 
         [0048]    The objective function value corresponding to the optimal solution indicates the minimal amount of delivery discrepancy using the MLC. Therefore, this value can be used to determine whether to pause the treatment beam to trade efficiency for higher dose conformality. By the same token, the optimization principle of the current invention can be naturally extended to incorporate other levels of tradeoff and result in a benign hybrid optimization problem. Additional physical considerations such as limited leaf velocity can be naturally incorporated into the invention by modifying the optimization setup with proper constraints. In principle, the underdose/overdose unitary cost needs to be assigned based on tissue property (tumor vs. normal tissue, parallel or serial, radiosensitivity, etc). In practice, such information is often known approximately. Despite this uncertainty, assigning inhomogeneous underdose and overdose unitary costs guides the sequencing algorithm towards a leaf configuration that offers improved tumor coverage and normal tissue protection. Accumulative dose discrepancy may also be incorporated by assigning a higher underdose unitary cost to voxels that have received lower dose than desired up to the time point of consideration, and conversely assigning a higher overdose unitary cost to voxels that are already overdosed. This mechanism would have a negative feedback feature to potentially prevent systematic underdose and/or overdose. 
         [0049]    The current invention includes adapting an optimized plan aperture to the instantaneous motion. This setup alleviates the burden of MLC modeling such as accounting for the tongue and grove effect, as they are implicitly incorporated in the treatment planning process. Similarly, starting from a treatment plan also relieves one of the task of sequencing MLC leaves and optimizing their velocities to achieve a desired intensity modulation, and allows one to focus on investigating the adaptation of an aperture under deformation. 
         [0050]    Two components contribute to the overall computation complexity of the proposed algorithm: (1) the numerical integration to yield c i  (or l i  in the homogeneous cost case) and (2) the search routine for the optimal placement. For the general inhomogeneous case, the numerical integration has a complexity proportional to the computation resolution in the BEV (Nx×Ny), where Nx and Ny are the discretization levels along and perpendicular to the leaf track respectively. The semi-exhaustive search scheme for optimal leaf placement has a complexity of ˜O(N x   2 ). With the relatively thin leaves and the overdose and underdose cost to be reasonably smooth, it is often the case that Ny&lt;&lt;Nx and the overall computation is dominated by ˜O(N x   2 ). Homogeneity of the unitary overdose/underdose cost and/or other structural assumptions (such as local convexity of the aperture), when utilized wisely, may further reduce the computation complexity. 
         [0051]    In summary, the current invention is a new optimization-based leaf adaptation methodology that explicitly penalizes underdose/overdose error. This method is generally applicable to all types of geometrical and topological changes of the anatomy, as well as various plan types. 
         [0052]    The present invention has now been described in accordance with several exemplary embodiments, which are intended to be illustrative in all aspects, rather than restrictive. Thus, the present invention is capable of many variations in detailed implementation, which may be derived from the description contained herein by a person of ordinary skill in the art. For example it is applicable when the initial plan is optimized with respect to a MLC geometry set that is different from the delivery MLC, (e.g., superresolution plan or conventional plan to be delivered with a finer MLC). Tongue and groove effect as well as delivery leakage can be naturally incorporated via substituting the indicator functions described above to more specific leaf geometry descriptors. 
         [0053]    All such variations are considered to be within the scope and spirit of the present invention as defined by the following claims and their legal equivalents.