Incorporating multiple targets in trajectory optimization for radiotherapy treatment planning

Abstract

A method for treatment trajectory optimization for multi-target radiotherapy treatment includes: determining a beam direction view (BEV) region and a BEV region connectivity manifold for each target group in a plurality of target groups, respectively. Information contained in the BEV regions and the BEV region connectivity manifolds for all target groups is used to guide an optimizer to find an optimal treatment trajectory. In order to improve the visibility of underexposed voxels of a planning target volume (PTV), a post-processing step may be performed to enlarge certain BEV regions that are considered for exposure during treatment trajectory optimization.

Description

[0001] This application is a divisional application filed on December 23, 2019, the international application date, entered the Chinese national phase on June 25, 2021, with application number 201980086475.8, and the invention name is “Incorporating multiple targets in trajectory optimization for radiotherapy treatment planning”. Background Art

[0002] Significant progress has been made in the inverse treatment planning of external beam radiotherapy using, for example, IMRT and VMAT treatment modes. With the increase in plan quality requirements and requirements for clinic patient throughput, the role of automated and more personalized treatment plans becomes increasingly important. Before the MLC leaf sequence and dose rate in the radiotherapy plan are inversely optimized, it may be necessary to select a treatment field, such as a VMAT trajectory or an IMRT field. The selection of treatment geometry (e.g., one or more isocenters, starting gantry angles and stopping gantry angles, and one or more collimator angles of the VMAT arc, or the aperture position, gantry angle, bed angle, and collimator angle of the IMRT field) may be determined by the clinical plan for a given treatment site. Considering the variability of the patient's anatomical structure and clinical goals, this plan may be suboptimal. In addition, radiotherapy treatment may involve tumors or lesions that may have very different volumes. Some tumors in the tumor may overlap partially or completely in space, while some of them may not intersect in space. Therefore, an improved method for optimizing treatment geometry is needed. Summary of the invention

[0003] In one aspect, the invention provides a trajectory optimization method for radiotherapy treatment of multiple targets.Optional features are defined in the dependent claims.

[0004] In another aspect, the invention provides a method for trajectory optimization of radiotherapy treatment of multiple targets.Optional features are defined in the dependent claims.

[0005] The present invention is a computer program and system for performing trajectory optimization of radiation therapy treatments.

[0006] According to some embodiments, methods for optimizing treatment geometry consider dose response to a region of interest (ROI). The ROI may include one or more planning target volumes (PTV) and one or more organs at risk (OAR). These methods can automate a large amount of manual planning work and can provide a new field of dose carving that can significantly enhance target dose consistency and reduce toxicity to normal tissues.

[0007] In some embodiments, beam direction view (BEV) regions and BEV region connectivity manifolds are determined by evaluating the dose response of each ROI at each vertex in the delivery coordinate space (DCS). The information contained in the BEV regions and BEV region connectivity manifolds is used to guide the optimizer in finding the optimal field geometry in the radiation therapy plan. According to some embodiments, in order to improve the visibility of underexposed voxels of the planning target volume (PTV), a post-processing step referred to herein as "coverage enhancement" can be performed to enlarge certain BEV regions that are considered for exposure during treatment geometry optimization.

[0008] In some embodiments, to optimize the treatment trajectory of a radiotherapy treatment of multiple targets, the BEV regions and BEV region connectivity manifolds are determined for each of the multiple target groups. The BEV regions associated with all target groups constitute a complete set of BEV regions, which are then used for treatment trajectory optimization. By determining the BEV regions for each target group separately, a more balanced plan for complex multi-target treatment can be achieved.

[0009] In some embodiments, the information contained in the BEV region and the BEV region connectivity manifold is used to guide the optimizer to find the best field geometry in the IMRT treatment plan. The treatment field can be represented as a set of nodes in the optimization framework. Beam angle optimization (BAO) may include: finding the set of nodes with the best maximum distance.

[0010] These and other embodiments of the present disclosure are described in detail below.For example, other embodiments relate to systems, devices, and computer-readable media associated with the methods described herein.

[0011] The nature and advantages of the embodiments of the present disclosure may be better understood with reference to the following detailed description and accompanying drawings. BRIEF DESCRIPTION OF THE DRAWINGS

[0012] The patent or application file contains at least one drawing in color. The International Bureau will provide copies of the patent or patent application publication with color drawings upon request and payment of the necessary fee.

[0013] Figure 1 is a schematic perspective view of a radiation therapy system.

[0014] Figure 2 is a schematic side view of a radiation therapy system.

[0015] Figure 3 A photon collimation system in a radiation therapy system is schematically shown.

[0016] Figure 4 An exemplary multi-leaf collimator (MLC) plane is shown.

[0017] Figure 5 Shows Figure 1 and Figure 2 Block diagram of an external beam radiation therapy system.

[0018] Figure 6 A representation of the delivery coordinate space (DCS) truncated to avoid collisions is shown, where each point (θ gantry ,θ couch ) is transformed into the physical three-dimensional (3D) position of the treatment head.

[0019] Figure 7 Examples of fiber bundle cross-section histograms of a differential beam direction view (BEV) and an integrated BEV fiber bundle cross-section histogram are shown in accordance with some embodiments.

[0020] Fig. 8A and Figure 8B An exemplary image showing BEV scores for two adjacent delivery coordinate vertices for chest wall treatment, respectively, according to some embodiments.

[0021] 9A to 9C Illustrated are exemplary steps for calculating a threshold BEV score for local variation, according to some embodiments.

[0022] FIG. 10A to FIG. 10B An example of the breast and supraclavicular node treatment areas (represented by the gold areas) before applying coverage enhancement according to some embodiments is shown ( Fig. 10A ) and an example of the area after applying coverage enhancement ( Fig. 10B ).

[0023] Fig.11 Illustrated is an exemplary quadrangularized spherical cube that may be used to evaluate angular flux in accordance with some embodiments.

[0024] Fig. 12A Illustrated are boundary nodes in a one-dimensional (1D) delivery coordinate space according to some embodiments.

[0025] Fig. 12B Illustrated are boundary nodes in a two-dimensional (2D) delivery coordinate space according to some embodiments.

[0026] Fig.13 A flow chart illustrating a method of trajectory optimization for radiation therapy treatment using slices with coverage enhancement in accordance with some embodiments is shown.

[0027] FIG. 14A to FIG. 14C Exemplary BEV maps of regions at the same vertex of three planning target volume (PTV) sets are respectively shown according to some embodiments.

[0028] Fig.15An overlay of four regions corresponding to three different PTV groups is shown in accordance with some embodiments.

[0029] Fig.16 A flow chart illustrating a method of trajectory optimization for radiation therapy treatment of multiple targets in accordance with some embodiments is shown.

[0030] Fig.17 A flow chart illustrating a method of toe angle optimization using slices according to some embodiments is shown.

[0031] Fig.18 A flow chart illustrating a method of toe angle optimization using slices according to some embodiments is shown.

[0032] Fig.19 A flow chart illustrating a method of beam angle optimization using slices with coverage enhancement in accordance with some embodiments is shown.

[0033] Fig. 20 A block diagram of an example computer system that may be used with systems and methods according to embodiments is shown. DETAILED DESCRIPTION

[0034] The present disclosure relates generally to treatment planning for radiotherapy using an external beam radiotherapy system, and more specifically to optimizing trajectories and field geometries in radiotherapy treatment plans. Beam direction view (BEV) regions and BEV region connectivity manifolds can be determined by evaluating the dose response of each region of interest at each vertex in a delivery coordinate space (DCS). The information contained in the BEV regions and BEV region connectivity manifolds can be used to generate optimized trajectories or optimized field geometries in radiotherapy treatment plans.

[0035] I. Treatment System

[0036] Figure 1 and Figure 2 Depicted is a radiation therapy system of the type that may be used in conjunction with the present invention. Figure 1 , a perspective view of a radiotherapy system (in this case a linear accelerator) is shown. Typically, such a system is capable of generating an electron (particle) beam or an X-ray (photon) beam for radiotherapy treatment of a patient on a treatment couch 35. Other radiotherapy systems are capable of generating heavy ion particles, such as protons. For the purposes of this discussion, only X-ray radiation is discussed. However, those skilled in the art will appreciate that the same principles apply to other systems.

[0037] The support 10 supports a rotatable gantry 20 with a treatment head 30. A control unit (not shown) is arranged next to the support 10, and the control unit includes a control circuit for controlling different operating modes of the accelerator. A high voltage source is provided in the support or in the gantry to provide voltage to an electron gun (not shown) on an accelerator guide located in the gantry 20. Electrons are emitted from the electron gun into a guide (not shown) in which they are accelerated. The source supplies RF (microwave) power to generate an electric field in the waveguide. The electrons emitted from the electron gun are accelerated in the waveguide by the electric field and leave the waveguide as a high energy (usually megavolt energy) electron beam. The electron beam then strikes a suitable metal target, thereby emitting high energy X-rays in the forward direction.

[0038] Now, refer to Figure 2 , shows a slightly more detailed side view of a radiation therapy system of the type that may be used in conjunction with the present invention. A patient P is shown lying on a treatment couch 35. X-rays formed in accordance with the above description are emitted from a target in a treatment head 30 in a divergent beam 104. Typically, the radiation beam is directed perpendicular to the target. Figure 2 The patient plane 116 of the page is positioned approximately one meter from the X-ray source or target, and the axis of the gantry 20 lies on the plane 116 so that the distance between the target and the isocenter 178 remains constant when the gantry 20 is rotated. The isocenter 178 is located at the intersection between the patient plane 116 and the central axis of the beam 122. The treatment volume to be irradiated is positioned about the isocenter 178.

[0039] Figure 3 A photon collimation system 300 is schematically shown having an upper aperture 310 (i.e., a Y1 aperture and a Y2 aperture; the Y1 aperture is omitted for clarity), a lower aperture 320 (i.e., an X1 aperture and an X2 aperture), and a multi-leaf collimator (MLC) 330. The field size in plane 340 at the isocenter 178 is indicated. The upper aperture 310, the lower aperture 320, and the leaves 332 of the MLC 330 include x-ray blocking material and are positioned in the head 30 to define the width of the x-ray beam at the patient plane. Typically, the apertures 310 and 320 are movable and, when fully opened, define a maximum beam of approximately 40 cm x 40 cm at the patient plane 116. The MLC 330 is located at the exit of the head 30 to further shape the x-ray beam. Since its introduction in 1990, the MLC has become a standard feature of most radiation therapy systems. Current examples of MLCs sold by the assignee of the present invention use up to 120 individually controllable leaves (typically thin tungsten blades) that can be moved in and out of the x-ray beam under control of the system software.

[0040] Figure 4An exemplary MLC plane is shown having a plurality of blades 332 arranged in opposing pairs and an aperture 415 created by movement of selected blades. Radiation passes through the aperture 415 and is shaped by it. Thus, the MLC can be used to collimate x-rays for conformal treatment of tumors from various angles ("3D conformal") and for intensity modulated radiation therapy ("IMRT"), whereby different radiation doses are delivered to different portions of the treatment area. The treatment volume (i.e., the irradiated volume along the path of the x-ray beam close to the isocenter 178) is defined by the apertures 310 and 320, the blade sequence of the MLC 330, and the collimator angle (i.e., the angle at which the MLC 330 is located in the head 30). Some external radiation therapy systems may include multiple MLC layers. The multiple MLC layers may be located in different planes and at different collimator angles.

[0041] Figure 5 Shows Figure 1 and Figure 2 Block diagram of an external beam radiotherapy system 500 of the present invention. The radiotherapy system 500 includes a beam source 510, a beam aperture 520, a gantry 530, and a bed 540. The beam source 510 is configured to generate a therapeutic radiation beam. The radiation beam may include x-rays, particles, etc. The beam aperture 520 includes an adjustable multi-leaf collimator (MLC) 522, which is used to spatially filter the radiation beam. The bed 540 is configured to support and position a patient. The bed 540 can have six degrees of freedom, namely, translation offsets X, Y, and Z, as well as rotation, pitch, and yaw.

[0042] A gantry 530 that orbits around a bed 540 houses a beam source 510 and a beam aperture 520. The beam source 510 is optionally configured to generate imaging radiation as well as therapy radiation. The radiation therapy system 500 may also include an image acquisition system 550 that includes one or more imaging detectors mounted to the gantry 530.

[0043] The radiation therapy system 500 also includes a control circuit 560 for controlling the operation of the beam source 510, the beam aperture 520, the gantry 530, the bed 540, and the image acquisition system 550. The control circuit 560 may include hardware, software, and memory for controlling the operation of these various components of the radiation therapy system 500. The control circuit 560 may include a fixed-purpose hard-wired platform, or may include a partially or fully programmable platform. The control circuit 560 is configured to perform one or more steps, actions, and other functions described herein. In some embodiments, the control circuit 560 may include a memory for receiving and storing a radiation therapy plan that defines one or more control points of a treatment field. The control circuit 560 may then send control signals to various components of the radiation therapy system 500, such as the beam source 510, the beam aperture 520, the gantry 530, and the bed 540, to execute the radiation therapy plan. In some embodiments, the control circuit 560 may include an optimization engine 562 configured to determine the radiation therapy plan. In some other embodiments, the control circuit 560 may not include an optimization engine. In those cases, the radiation treatment plan may be determined by an optimization engine in a separate computer system, and then transmitted to the control circuitry 560 of the radiation treatment system 500 for execution.

[0044] II. Radiation therapy planning

[0045] Radiation therapy is typically implemented according to a radiation therapy plan, which typically takes into account the desired radiation dose specified to be delivered to the tumor and the maximum dose that can be delivered to surrounding tissue. Various techniques for developing radiation therapy plans can be used. Preferably, a computer system used to develop radiation therapy plans provides outputs that can be used to control the radiation therapy system (including control points and MLC leaf movements). Typically, the desired dose specified in the radiation therapy plan is delivered over several sessions (called fractions).

[0046] Several techniques have been developed to generate radiation treatment plans for IMRT or conformal radiation therapy. Typically, these techniques involve solving the "inverse" problem of determining the optimal combination of angles, radiation doses, and MLC leaf movements to deliver the desired total radiation dose to the target while minimizing irradiation of healthy tissue. The inverse problem is even more complex for planning arc therapies such as volumetric modulated arc therapy (VMAT), where one or more external treatment coordinates such as isocenter position, gantry angle, bed angle, and bed offset are in motion while the target volume is irradiated. To date, radiation oncologists or other medical professionals such as medical physicists and dosimetry physicians have used one of the available algorithms to develop and optimize radiation treatment plans.

[0047] Typically, such planning begins with volumetric information about the target tumor and any nearby tissue structures. For example, such information may include a map of a planning target volume ("PTV"), such as a prostate tumor, which is prescribed by a physician to receive a certain therapeutic radiation dose with an allowable tolerance. Volumetric information about nearby tissues may include, for example, maps of the patient's bladder, spinal cord, and rectum, each of which can be considered an organ at risk (OAR) that can only receive a much lower maximum prescribed radiation dose. This volumetric information, along with prescribed dose limits and similar goals set by medical professionals, is the basis for calculating an optimized dose distribution (also referred to as a fluence map), which in turn is the basis for determining a radiotherapy plan. For example, the volumetric information can be simplified to an objective function or a single quality factor to illustrate the relative importance of various trade-offs inherent in radiotherapy planning and the constraints that must be met to make the radiotherapy plan medically acceptable or physically possible.

[0048] III. Beam Emission View (BEV) Slices

[0049] A. Delivery Coordinate Space (DCS)

[0050] The most advanced technology for optimizing treatment trajectories (e.g., VMAT trajectories) in external beam radiotherapy treatments involves dosimetric characterization of candidate injection directions. Considering the dose response of both the PTV and OARs within the patient, the goal of optimization may be to identify which injection directions in the admissible delivery coordinate space are more suitable for treating the patient.

[0051] The delivery coordinate space (DCS) is the set of all permissible coordinates that parameterize the configuration of the delivery device and is truncated to avoid collisions (e.g., machine-to-machine and machine-to-patient). For a C-arm linear accelerator with a fixed isocenter, a point in the delivery coordinate space can be defined as (θ gantry ,θ couch ) tuple. A DCS can be discretized into a 2D grid defined by a collection of vertices (racks, bed values), edges, and triangular faces. Therefore, It can be represented by a simplicial complex (grid), which is defined as follows:

[0052] Vertex: (θ gantry ,θ couch ) angle tuple;

[0053] Edges: ordered pairs of start and end points; and

[0054] Triangle: An ordered list of 3 vertex indices.

[0055] Figure 6 A representation of the DCS is shown as a 3D grid, where each point 610 (θ gantry ,θ couch ) is transformed to the physical three-dimensional (3D) position of the treatment head, and the space is truncated to avoid collisions.

[0056] B.BEV Slicing

[0057] Reference [1] discusses a method for optimizing trajectories in radiotherapy using slices (referred to herein as TORUS). Dosimetry experience in radiotherapy planning has shown that BEV provides a valuable tool in determining the geometric setup for both dynamic gantry treatments (e.g., VMAT) and static gantry treatments (e.g., IMRT).

[0058] For the source position corresponding to the vertex v in the 3D DCS The BEV plane (also called the isocenter plane) can be defined as being perpendicular to the vector And including the isocenter The BEV plane can be discretized into a plane with N x ×N y A 2D array of pixels, where each pixel on the 2D grid represents a single beamlet.

[0059] To detect all possible beamlets, the intensity of each beamlet can be set to unity, and the 3D dose response of each beamlet can be evaluated. The 3D dose can be processed to determine the dose statistics for each region of interest (ROI) for each beamlet. ROIs can include, for example, planning target volumes (PTVs) and organs at risk (OARs). If the index n at position (x, y, z) is ROI ROI, delivering coordinates of vertex n v The pixels in a given BEV plane (n x , n y ) is given by If given, the index is n ROI The volume integrated dose of the ROI can be estimated as follows:

[0060]

[0061] If there is N ROI Region of interest, then the BEV dose bundle cross section The size is A 4D array containing the volume-integrated dose values ​​for each ROI from each beamlet.

[0062] In some embodiments, the BEV fractional cross section can be defined as a cross section from the 4D BEV dose bundle To size A reduction of a 3D matrix of , where the values ​​are a measure of the "goodness" of the beamlets based on the ROI dosimetry. In some embodiments, the "goodness" score can be evaluated as a linear combination of the dose for each ROI,

[0063]

[0064] coefficient It can be set by the user. For example, for the body, the value of the coefficient can be set to -0.2; for OAR, the value of the coefficient can be set to -1 to -10 (for example, more negative weights can be given to critical organs); and for PTV, the value of the coefficient can be set between 0 and +1.

[0065] C. BEV regional connectivity manifold

[0066] The BEV regional connectivity manifold can be constructed in two steps. First, the information contained in the BEV fractional bundle cross-sections can be considered and a binary selection process can be applied to determine whether a given pixel (sub-beam) is a "good" or "bad" treatment candidate. For each BEV plane, the set of "good" sub-beams forms a region. Each region consists of a set of continuous pixels in the BEV plane and represents a potential opening candidate for optimization. Next, it can be determined how the region is connected to other regions in adjacent vertices. The final structure consisting of regions and their connections forms the BEV regional connectivity manifold.

[0067] 1. BEV Regional Score

[0068] In some embodiments, if a beamlet intersects the PTV and its score S is above a certain threshold Then it can be considered a "good" candidate. Selecting the appropriate threshold can be a nontrivial task and may vary from case to case. For example, the threshold for a beamlet treating a superficial target with little body or OAR obstruction (e.g., prone breast irradiation) may be different from the threshold for a deep target (e.g., in prostate treatment) where the best possible approach may still be to penetrate healthy tissue to substantial depth for treatment.

[0069] According to some embodiments, the region score Can be used to define regions where potentially “good” beamlets are defined as having

[0070]

[0071] Where N PTVis a set of ROI indices for the PTV region of interest. Therefore, if a beamlet intersects the PTV and its score is above a certain threshold Then it can be considered a "good" candidate. The region score R can classify the sub-bundle into multiple regions, and can also be used as a normalized score of the goodness of the sub-bundle (eg, the maximum region score is 1).

[0072] 2. Determination of score threshold

[0073] According to some embodiments, the score threshold The histogram of the BEV fiber bundle cross section can be used to automatically determine the dose volume histogram (DVH) in the spirit of the dose-volume histogram. Where n is a set and the subset of indices considered The associated BEV fiber bundle cross-section histogram can be defined as follows:

[0074] Determine the set Maximum and minimum values ​​of: F max and F min .

[0075] Create a range from F min to F max N bins , and each N bins Initialized to 0. These can be called difference histogram bins, and for each bin index N bins , expressed as

[0076] ·Loop through each And make the value The corresponding bin is incremented by 1.

[0077] The integral histogram bins can be defined as

[0078] Normalize the difference histogram bins and the integral histogram bins to have a maximum value of 1.

[0079] Figure 7 An example of a differential BEV fiber bundle cross-section histogram 710 (dashed line) and an integrated BEV fiber bundle cross-section histogram 720 (solid line) according to some embodiments is shown. Dashed line 730 intersects the integrated histogram curve 720, and the vertical height to horizontal length ratio may be 0.7. Their intersection with the abscissa may give a score threshold

[0080] Using this definition, the BEV fiber bundle cross section can be (BEV fiber bundle cross-sections only represent the total PTV dose) Define BEV PTV dose histograms. These histograms can be expressed as and In some embodiments, the histogram can be used to determine a temporary PTV dose threshold so that only the fraction of beamlets with the top 50% highest PTV doses are considered. The threshold can be expressed as And for the integral histogram bins with height 0.5 dosage.

[0081] Next, the thresholded BEV score histogram can be used to limit the index BEV fraction bundle cross section These threshold histograms can be expressed as and The ratio of the vertical height to the horizontal height of a point on the integral histogram curve (for example, Figure 7 The side length ratio of the rectangle formed by the axis shown in and the dashed line 730 is as follows:

[0082]

[0083] in represents the integral histogram height at the score value S (i.e., where N bins is the corresponding bin index). The ratio is from S = S min +∞ to S=S max 0 (where B S =0) changes monotonically. Score threshold can be defined such that the ratio ratio(S threshold )=0.7.

[0084] It should be understood that the score threshold determination method described above is only an example. According to other embodiments, other determination methods may be used.

[0085] Fig. 8A and Figure 8B An exemplary image of the BEV scores for two adjacent delivery coordinate vertices for chest wall treatment is shown. Darker shading represents lower score values, while lighter shading represents higher score values. Gold represents beamlets that pass the region selection criteria and are therefore considered BEV regions. In this example, from the BEV perspective, Fig. 8A There are two disconnected BEV regions 810a and 820a in the BEV shown, and Figure 8B There are two disconnected BEV regions 810b and 820b in the BEV shown.

[0086] 3. BEV Regional Connectivity Manifold

[0087] Regions at adjacent vertices in the delivery coordinate space can be stitched together to form a complete BEV region connectivity manifold. The BEV region connectivity manifold contains information about how candidate target regions in the BEV change, appear, split, and disappear as one moves along the delivery coordinate space in all directions. For example, Fig. 8A and Figure 8B In the example shown, two images of two adjacent delivery coordinate vertices are overlapped, and regions 810a and 810b on the left can overlap each other, while regions 820a and 820b on the right can overlap each other. Therefore, it can be inferred that regions 810a and 810b are connected, and regions 820a and 820b are connected when moving along the edge in the delivery coordinate space. The set of all regions and all connections along the edge of the delivery coordinate space forms a BEV region connectivity manifold.

[0088] The information contained in the BEV regions and the BEV region connectivity manifold can be used to generate optimized trajectories or field geometries in radiotherapy. For example, as discussed below and in reference [1], trajectory optimization methods can use the BEV region connectivity manifold as a scaffold to guide the optimizer, which can make the search space small enough that graph search techniques with efficient computational time can be applied.

[0089] D. Coverage Enhancement

[0090] As discussed above, beamlets that intersect the PTV and have a BEV score above a global threshold BEV score may be deemed desirable and included in the BEV region, and these beamlets may be considered for exposure during treatment trajectory optimization. This process may not guarantee that every element volume (e.g., voxel) of the PTV can be irradiated from a sufficient number of incident directions.

[0091] According to some embodiments, to improve visibility of underexposed PTV voxels of the PTV, a post-processing step referred to herein as "coverage enhancement" may be performed to enlarge certain BEV regions. Coverage enhancement may be achieved by establishing coupling between voxels of each PTV group and BEV regions, as described below.

[0092] To preferentially enhance PTV coverage of BEV regions, a beamlet-to-beamlet variation of the threshold BEV fraction may be introduced by considering the 3D incident direction on each PTV. According to some embodiments, the threshold BEV fraction for local variation may be calculated as follows. Assume that the total number of beams under consideration is N b , then a coverage fraction p (e.g., 0.3) can be defined so that the minimum number of beams that can expose each PTV voxel is pN bEach PTV voxel can be projected onto the BEV score map for each bundle under consideration, and its score distribution can be calculated. b So the exposure pN can be determined using the voxel-wise BEV score distribution. b The voxel-specific threshold BEV score for each bundle is th voxel, the voxel-specific threshold BEV score is denoted as S′ i Assume S threshold is the original global threshold, then the true voxel-specific threshold is defined as S i =min(S′ i , S threshold ).

[0093] These voxel-based score thresholds S can then be used i To define a threshold BEV score for local variation. First, voxels are projected as points onto each BEV, and pixels inside their mosaic are then defined using linear interpolation, and pixels outside the mosaic are defined using nearest neighbor interpolation. To smooth artifacts introduced by pixel grid inconsistencies, a small amount of Gaussian smoothing can be applied to obtain a final score threshold map for local variation. Regions can then be defined as sets of pixels that intersect the PTV and have a score above their local threshold.

[0094] 9A to 9C Exemplary steps in the calculation of the score threshold of local variation as described above are illustrated. Fig. 9A The image in shows the result of applying linear interpolation inside the PTV voxel mosaic (0.0 corresponds to the original global threshold); Fig. 9B The image in shows the result including outer nearest neighbor interpolation and Gaussian smoothing; and Fig. 9C The resulting increase in region shape when a new local variation score threshold is used is illustrated. Fig. 9C The blue area 920 in FIG. 9 indicates the BEV area before coverage enhancement, while Fig. 9C The pink area 920a in FIG. 8 represents the enlarged BEV area after coverage enhancement.

[0095] FIG. 10A to FIG. 10B An example of a BEV area 1010 (indicated by the gold area) during breast and supraclavicular node treatment before coverage enhancement is applied is shown ( Fig. 10A ) and an example of region 1020 after applying coverage enhancement ( Fig. 10B As shown, the BEV area 1020 after coverage enhancement is enlarged compared to the BEV area 1010 before coverage enhancement. For example, (in this case, corresponding to the smaller supraclavicular node) Fig. 10BThe sub-region 1020a shown is not included in the BEV region 1210, so those voxels may not be "hit" by the beamlet. After coverage enhancement, these voxels may be fully exposed to radiation.

[0096] IV. Trajectory Optimization for Radiotherapy Treatment

[0097] Reference [1] discusses trajectory optimization in radiotherapy using slices (referred to herein as TORUS). The TORUS method uses BEV regions and BEV region connectivity manifolds as guidance to automatically generate heuristic optimal radiotherapy trajectories to effectively deliver high-quality VMAT treatment plans. TORUS uses a dual-metric optimization to generate optimal treatment trajectories using an optimization graph on the delivery coordinate space. The nodes in the optimization graph can represent individual control points, and the trajectory can be defined as the path that minimizes the minimum distance metric, while the maximum distance metric can serve as a measure of the goodness of selecting the optimal trajectory.

[0098] A.PTV angular flux

[0099] One of the concepts used in the TORUS approach is PTV angular flux, which is related to the novelty of the three-dimensional (3D) directional vectors of the incident beamlets at a given point in the PTV. Inverse dose optimization may perform better at more angles at which the radiation beam enters the patient. The reason for this may be twofold. First, by entering the patient from multiple directions, the ratio of overlapping dose within the PTV to the surrounding OARs may be greater, resulting in a steeper dose gradient outside the PTV. Second, each beamlet from each direction may provide a different 3D dose contribution to the patient. Therefore, increasing the number of such unique beamlets may give the optimizer more "basis vectors" to work with when carving an optimal dose distribution around critical structures.

[0100] Note that simply entering the patient from multiple directions may not be sufficient to ensure optimal plan quality. There may be some cases where, even with a large number of beams entering the patient, portions of the PTV may only be exposed from a few directions while protecting nearby OARs. This may result in insufficient coverage of small areas of the PTV, non-conformal areas (dose streaks), or poor dose compensation. To encourage maximum coverage and conformality, it may be desirable to target each element volume of the PTV individually from many different directions.

[0101] According to some embodiments, the angular flux at a given point in the PTV may be estimated by computing the 3D direction vectors of the incident beamlets and binning them in angular bins. Fig.11An exemplary quadrangularized spherical cube that can be used to evaluate angular flux is illustrated. As illustrated, each cube face can be divided into 4×4 squares. Thus, there are 6×4×4=96 squares. Each square corresponds to a single angular bin. This can provide a bin size of approximately 20 degrees, which can correspond to the same order of magnitude distance between vertices in the delivery coordinate space (DCS). In general, each cube face can be divided into 2 2n squares, where n is a positive integer. Fig.11 In the illustrated example, n=2. This binning approach may result in bins of unequal solid angles, but the difference may be relatively small (up to 19%). Additionally, randomly orienting each spherical cube at a set of sampling points may offset this difference. In other embodiments, the angular flux may be evaluated by calculating the 3D direction vectors of incident beamlets passing through a closed surface of a different shape than the quadrangularized spherical cube. For example, the closed surface may be a spherical surface or a cubic surface.

[0102] The set of sampling points may be distributed within the PTV. During trajectory optimization, the angular flux of each sampling point may be evaluated and optimized. According to some embodiments, the angular flux of a given PTV point may be stored as a 12-byte bit set, indexing each bin from 0 to 95, thereby enabling fast bit-by-bit calculations to be performed. If n PTV points are considered, the angular flux state is represented as a vector of bit sets of length n, each bit set having 96 bits.

[0103] B. Dual distance metrics

[0104] According to some embodiments, the information stored in the BEV fractional bundle cross-sections and the BEV regional connectivity manifold can be used to generate a treatment plan for radiotherapy. The BEV regional connectivity manifold can be used as a scaffold for guiding the optimizer, which can make the search space small enough to apply graph search techniques with fast computation time.

[0105] In general, the optimal trajectory in VMAT-type treatments may be one that hits the PTV as much as possible, avoids or minimizes healthy tissue dose, and enters the PTV from many different directions, and can be accomplished in a relatively short delivery time. Some of the goals may conflict with each other. For example, a treatment plan that treats each PVT element from every direction with a high degree of MLC modulation may produce a nearly ideal dose distribution, but may also take too much time to deliver. Therefore, it may be desirable to find a trajectory that covers many directions around the patient in an efficient manner. According to some embodiments, an optimization method may seek to maximize the trajectory length while using a relatively "straight" trajectory for efficient delivery. The "straightness" of a trajectory can be understood in terms of geodesics. A geodesic is the shortest path between two given points in a curved space. A geodesic can be calculated by finding the line that minimizes the distance function between the two points.

[0106] According to some embodiments, based on the BEV fractional bundle cross-section and the BEV regional connectivity manifold, there are regions that traverse high regional scores. Trajectories with control points of may be preferred (e.g., as defined by equation (3)). To encourage geodesics to traverse these control points, a distance function can be defined such that small non-negative distances are preferred. On the other hand, to select long trajectories, another distance function can be defined such that larger distances are preferred.

[0107] To overcome this inherent conflict between trying to find a "short" trajectory and a "long" trajectory, according to some embodiments, a dual distance metric scheme is used in the optimization. The dual distance metric includes two distance functions that play different roles in the optimization. These two distance functions can be referred to as the minimum distance function D min (where smaller values ​​are preferred; minimization defines the geodesic path in the graph) and the maximum distance function D min (where larger values ​​are preferred; maximization is achieved by selecting trajectories).

[0108] C. State diagram optimization

[0109] According to some embodiments, a symmetric directed graph can be used to use a minimum distance function D min and the maximum distance function D max Perform trajectory optimization. Since the distance function may depend not only on the edges of the graph, but also on the history of the trajectory to reach that point, statefulness can be introduced into the graph.

[0110] The graph can include a set of nodes and a set of edges connecting these nodes. To take into account the file history in the graph, a vertex can be defined as a node and state pair. This state can be the PTV angular flux state described above. And edge E = (N1, N2), the successor point is The new status The successor state function Given, the successor state function describes how the state moves along the edge E. The minimum and maximum distances between these two points can be expressed as D min / max (P1, P2).

[0111] A path can be an ordered sequence of points P(P1, ..., P n ), so that for all 1≤i<n, and The minimum and maximum distances for this path may be Given an initial state and two nodes N and N′, assuming that the set of all possible paths between them is P (P1, ..., P n ), so that and in Not specified, represented by The trajectory set can be defined as the set of paths with the minimum distance,

[0112]

[0113] Then, the optimal set of trajectories between these two nodes can be defined as the trajectory with the largest distance,

[0114]

[0115] Using the definitions provided by equations (5) and (6), the minimum and maximum distances between two arbitrary nodes are in

[0116] Given an initial state Now the goal can be to define the global optimal trajectory. If attention is limited to a certain set The start and end nodes in the set The optimal set of trajectories ending on can be defined as follows:

[0117]

[0118] in in

[0119] The optimization problem may be for an initial state and a collection of start and end points Find the element

[0120] D. Graphics Definition

[0121] According to some embodiments, a control point in the graph can be uniquely determined by three integers: a vertex v (i.e., a point on the delivery coordinate space or DCS), a collimator index c (which determines the collimator angle from a discrete set of possibilities), and a region bitfield b. The region bitfield b is a list of Boolean flags that determine which region subset is selected for a given vertex. These three integers can define a node in the graph as N = (v, c, b). The starting MLC blade position can be determined by fitting to this region subset of the BEV.

[0122] According to some embodiments, the BEV region connectivity manifold may include multiple mutually disjoint connected components, each of which forms a single connected portion of the total search graph by the following edge definition. Given two nodes N1 = (v1, c1, b1) and N2 = (v2, c2, b2), there may be an edge E connecting these two nodes if the following conditions are met:

[0123] There is an edge connecting v1 and v2 in the delivery coordinate space;

[0124] .Δt collimator ≤Δt directional , where Δt collimator is the collimator movement Δθ collimator =θ collimator (c2)-θ collimator (c1) time, and Δt directional is the time to move from vertex v1 to v2 in rack and bed space.

[0125] Boundary node set are potential starting nodes and auxiliary nodes for graph optimization. To define this set of boundary nodes, it may be necessary to first define boundary vertices in the one-dimensional (ID) space and two-dimensional (2D) space (which can be generalized to higher dimensions) in the delivery coordinate space. The 1D space may consist of only vertices and edges, and the boundary vertices are vertices that touch at most a single edge. Similarly, in the 2D space of vertices, edges, and faces, the boundary vertices are vertices that belong to edges that touch only a single face. According to this definition, the region is Possibly a border region:

[0126] (1) Index n v The vertex v is the boundary vertex;

[0127] (2) There is a contact index n in the delivery coordinate space vThe edge of vertex v is such that there is no region edge emanating from region r along this delivery coordinate space edge.

[0128] FIG. 12A to FIG. 12B These conditions are illustrated graphically. Fig. 12A A 1D delivery coordinate space may be represented where boundary vertices 1210 and 1220 are represented as solid circles. Fig. 12B A 2D delivery coordinate space can be represented. The BEV region connectivity manifold is represented by ellipses (each of which represents a region) in the five BEV planes 1230a to 1230e and the connections between them. Region 1240 in the leftmost plane 1230a and region 1250 on the rightmost plane 1230e are boundary regions that satisfy the above condition (1), while region 1260 in plane 1230d has no connection to the region on the right and is therefore a boundary region that satisfies the above condition (2). According to these definitions, a node is a boundary node if all regions represented by its region bit field are boundary regions.

[0129] The state information used in graphics optimization is PTV angular flux Using the bit set definition of the PTV angular flux as described above, the successor function σ defining how the state changes can be obtained by the bitwise OR operator Definition, where is the contribution to the diagonal flux state from the region represented by the region bitfield b2 at vertex v2.

[0130] Given two points Where N i =(v i , c i , b i ), i = 1, 2, edge E connects nodes, and the minimum distance function / maximum distance function can be defined as follows:

[0131]

[0132] D max (P1, P2)=Δθ(E)S(P1,P2), (9)

[0133] Where Δθ is the physical angular distance traveled by the treatment head. If the resulting MLC configuration violates machine limitations, the minimum distance function is set to +∞. Note that the fractional terms S(P1, P2) appear as reciprocals of each other in each equation, reflecting the fact that, roughly speaking, maximum distance = "goodness" and minimum distance = 1 / "goodness". The definitions and meanings behind each term are as follows:

[0134] ·

[0135] - in is the average area fraction of bit field b regions per unit area, A(b) is the combined area of ​​bit field b regions, and A penalty (b) Total non-region area exposed by the fitted MLC. This term penalizes poor MLC target fits and encourages high scoring regions.

[0136] - is a term that encourages regions that contribute novel angles to the existing angular flux state. If is the average contribution (number of bits) of each region to the blank angular flux state, then for or 1 (whichever is greater), which is divided by Normalization.

[0137] ·t mlc is the time it takes the MLC to move between control points. This penalizes collimator angles that cause excessive MLC motion.

[0138] ·p eff =(1-Δθ / Δθ max ) is a factor that penalizes the edges of the trajectory where they are almost stationary, so that no treatment time is wasted at these locations.

[0139] The fractional term S(P1, P2) is likely the main force driving the optimizer to find a trajectory that targets the PTV from a good direction and gives contributions from different directions while avoiding poor MLC target fits.

[0140] E. Graphics Optimization Solution

[0141] The state dual distance metric graph optimization defined above can be solved using Dijkstra's algorithm. Converting the list of points from the resulting optimal trajectory to control points can provide the desired radiation therapy trajectory. Running Dijkstra's algorithm from the start node without stopping at any particular end node may produce a tree structure of points, effectively completing a 1 to N search from a given start node to all other nodes, where the computational complexity is the same as a normal 1 to 1 search between two nodes. By selecting the trajectory with the maximum distance as the largest, the restricted optimization problem of finding the best trajectory from a given node can be effectively solved.

[0142] In general, to find the global best trajectory with the largest maximum distance, it may be necessary to repeat this calculation from each possible node (N to N search). An approximately optimal trajectory can be found by using the end node of the previous run as the start node of the consecutive run, picking an arbitrary start node and repeatedly running the algorithm with the same initial state. In some embodiments, this process is repeated twice; therefore, trajectory optimization can be performed with the same computational complexity as the underlying Dijkstra algorithm.

[0143] For each connected component of the graph, a path optimization can be performed. Finally, the trajectory with the largest maximum distance among all possibilities can be selected. The presence of angular flux states in the distance function can help ensure that the selected trajectory will also be one that tends to provide novel directions from which the PTV is treated.

[0144] F. Trajectory Optimization Methods for Radiotherapy Treatment Using Slices with Coverage Enhancement

[0145] Fig.13 Shown is a flow chart illustrating a method 1300 of trajectory optimization for radiation therapy treatment using slices with coverage enhancement in accordance with some embodiments.

[0146] At 1302, a patient model is provided. The patient model includes one or more planning target volumes (PTVs), the one or more PTVs including a first PTV.

[0147] At 1304, a delivery coordinate space (DCS) is defined. The DCS may include a plurality of vertices. Each respective vertex defines a respective beam direction view (BEV) plane.

[0148] At 1306, for each respective BEV plane of the respective vertex and for each respective pixel of the respective BEV plane, a dose from a first PTV of a respective beamlet originating from the respective vertex and passing through the respective pixel of the respective BEV plane is evaluated. A respective BEV score for the respective pixel of the respective BEV plane is also evaluated based at least in part on the dose from the first PTV of the respective beamlet.

[0149] At 1308, a global threshold BEV score for the first PTV is determined based on the BEV scores of the pixels of the BEV plane.

[0150] At 1310, for each respective BEV plane of the respective vertex, one or more initial BEV regions are determined by comparing the respective BEV score of each respective pixel to a global threshold BEV score. The BEV scores of the pixels within the one or more initial BEV regions are greater than or equal to the global threshold BEV score. Each respective initial BEV region may expose one or more voxels of the first PTV to a beamlet corresponding to the respective vertex.

[0151] At 1312 , for each respective voxel of the first PTV, a voxel-specific threshold BEV fraction is determined based on the initial BEV area of ​​the BEV plane such that the respective voxel is exposed by beamlets corresponding to a predetermined proportion of vertices out of the total number of vertices.

[0152] At 1314, for each respective BEV plane of the respective vertex, and for each respective pixel of the respective BEV plane, a local threshold BEV score for the respective pixel is determined based on the voxel-specific threshold BEV score. One or more BEV regions of the respective BEV plane are determined by comparing the respective BEV score of each respective pixel to the local threshold BEV score. The BEV scores of the pixels within the one or more BEV regions are greater than or equal to the local threshold BEV score.

[0153] At 1316, a BEV region connectivity manifold of the first PTV is determined. The BEV region connectivity manifold represents connections between BEV regions of adjacent vertices.

[0154] At 1318, one or more optimal treatment trajectories are selected based on the BEV region connectivity manifold of the first PTV.

[0155] You should understand that Fig.13 The specific steps shown provide a specific method for trajectory optimization for radiotherapy treatment according to some embodiments. Other step sequences may also be performed according to alternative embodiments. For example, alternative embodiments may perform the steps in different orders. Moreover, each step may include multiple sub-steps, which may be performed in various orders depending on the circumstances of the individual steps. Further, additional steps may be added and some steps may be removed depending on the specific application. Those of ordinary skill in the art may recognize many variations, modifications, and alternatives.

[0156] V. Trajectory Optimization for Multi-Target Radiotherapy Treatment

[0157] Radiation therapy treatments may treat multiple tumors or lesions that may have very different volumes. Some of these tumors may partially or completely overlap spatially, while some of them may not be spatially intersecting. An example of full spatial overlap may be enhancement targets within the primary planning target volume (PTV). For example, partial overlap may occur in head treatment and neck treatment. Spatially non-intersecting lesions may be encountered as multiple intracranial metastases. For example, breast cancer treatment may involve multiple targets defined for breast tissue and associated nodal groups.

[0158] The original TORUS scheme treats all PTV volumes as a single entity, where target-specific features are not taken into account. For example, no record of enhanced regions within the PTV is kept. In addition, targets with small volumes are not considered to be on an equal footing with targets with large volumes. According to some embodiments, the TORUS scheme can be extended for radiotherapy treatment of multiple targets. In the following, various aspects of multi-target trajectory optimization are discussed.

[0159] A. Sections of specific regions of the PTV group

[0160] According to some embodiments, in a multi-target radiotherapy treatment, the region of each target may be calculated independently while setting the weights of other PTVs to zero. By calculating the score threshold for each target separately, a more balanced approach to complex multi-target treatments may be achieved. During optimization, the region associated with each target of each BEV may constitute the entire set of regions. To ensure that each target is treated fairly in the optimization, the integrated region score associated with each target and the number of points sampled in the target to calculate the PTV angular flux novelty may be proportional to an overall weight factor that is associated with the target.

[0161] In some embodiments, the identity of each PTV group can be kept intact relative to the identities of other PTV groups. Without loss of generality, it can be assumed that N PTV The target structures are divided into N groups, where 1≤N≤N PTV In some embodiments, the target set may correspond to a volumetric region to be associated with an optimization objective. In some other embodiments, the target set may include a plurality of partially spatially overlapping volumes or fully spatially overlapping volumes.

[0162] BEV Area A BEV map at each vertex of the Discretized Delivery Coordinate Space (DCS) for each PTV group indexed i=1, ..., N may be included. A vertex is a tuple where r iso is the isocenter position, and s a are generalized coordinates corresponding to the mechanical axes of the treatment device. For example, in the case of a C-arm linear accelerator and a constant isocenter, the vertex may take the form v C-arm =(r iso ,θ gantry ,θ couch ). As discussed above, the BEV region It may be indicated which areas of the BEV map through which PTV group i is to be irradiated. There may be one BEV map for each PTV group at each vertex.

[0163] FIG. 14A to FIG. 14C The BEV area at the same vertex of three PTV groups (i=1, 2, 3) is shown. . Gold areas 1410, 1420, 1422, and 1430 indicate BEV areas through which radiation can be administered to a target. Information about the BEV areas 1410, 1420, 1422, and 1430 of each BEV may be stored in the associated area Dark regions 1414, 1424, and 1432 may correspond to regions through which irradiation may deposit too little dose to the target group or too much expense to healthy tissue or organs at risk (OARs). In some embodiments, the region fractions of BEV regions 1410, 1420, 1422, and 1430 may be in the range (0, 1], while the region fractions of dark regions 1414, 1424, and 1432 may be in the range (-∞, 0] (e.g., as expressed in equation (3) above).

[0164] exist FIG. 14A to FIG. 14C In the example shown, Fig.14A The BEV area 1410 in the BEV map shown may correspond to the primary PTV (i=1), and Fig. 14B and Fig. 14C BEV regions 1420, 1430, and 1440 in the BEV map shown in (i=2, 3) may correspond to sub-PTVs within the main PTV to which two different magnitudes of boost doses may be delivered. These PTV groups are examples of complete spatial overlap between PTVs; that is, Fig. 14B and Fig. 14C The enhancement volume shown is located in volume Fig.14A The interior of the main PTV is shown.

[0165] BEV areas of all vertices of all PTV groups The entire set of regions may be constructed. When generating trajectories, the generalized TORUS algorithm may decide which BEV regions 1410, 1420, 1422, and 1430 of which target groups are to be included in the MLC aperture at each incident direction.

[0166] B. PTV group specific corner density

[0167] As described above, one of the concepts in the TORUS method is the PTV angular flux, which characterizes how a set of element volumes (i.e., voxels) of a PTV are irradiated via a given beamlet. In the presence of multiple PTV groups, the generalized TORUS method may strive to irradiate each element volume of each PTV in each PTV group from as many directions as possible.

[0168] In some embodiments, the weighted volume of each PTV group It can be calculated as follows:

[0169] v j (x) = n j (x)ΔV(x); (10)

[0170]

[0171]

[0172] where n j (x) is the relative volume of voxel x occupied by PTV j, ΔV(x) is the volume of voxel x, and the set contains the index of the PTV in group i, and ω j (x) is a non-negative PTV-specific weight. Then, each PTV group i can be assigned a corner density n A，i , the corner point density n A，i Written as:

[0173] N A (V) = αV β +λ; (13)

[0174]

[0175] Where N A (V) is a function that returns the number of sampling points, and α, β, and λ are parameters that can be tuned to maintain a sufficiently large point density over all relevant volumes V. In some embodiments, β can be a positive number less than one.

[0176] In some embodiments, the PTV-specific weight ω j (x) can be a positive fractional coefficient, for example, These positive fractional coefficients are used to calculate BEV pixel specific fractional region values. In some other embodiments, different relative weightings between PTV groups may be employed.

[0177] C. PTV group-specific region connectivity manifold

[0178] According to some embodiments, based on A set of regional connectivity manifolds can be constructed Keep track of how and which regions in PTV group i in different connected vertices are connected to each other. If the configuration of the delivery device allows, there may be edges (ie, connections) between neighbor vertices. Maintain an ordered list of all regions through which PTV group i is to be irradiated And the vertex set of DCS A list of specific vertices in the region in

[0179] D. Collapse of PTV-specific BEV maps into a single multi-PTVBEV map

[0180] As discussed above and in reference [1], TORUS uses a graph optimization approach that is a modification of Dijkstra's algorithm that outputs an entire tree of minimum distances from a starting node to all other nodes on the graph. According to some embodiments, in order to reduce the computational workload in the case of multiple PTV groups, a strategy may be adopted to reduce or collapse the PTV group-specific BEV maps onto a single multi-PTV BEV map. In one embodiment, the collapse may be performed on a BEV pixel-specific basis by assigning the pixel value to the multi-PTV BEV map at the pixel coordinate (x, y) of vertex v,

[0181]

[0182] In other embodiments, the collapse may be performed using other processes that preserve the property that if for some i, Then it remains within (0, 1] on the collapsed map.

[0183] Fig.15 An overlap of four regions 1510, 1520, 1530, and 1540 corresponding to three different PTV groups is shown. Regions 1520 and 1530 shown in purple and region 1540 shown in green correspond to sub-regions in which a collapse function has been used to assign region values ​​to affected pixels in the BEV. Due to the collapse, the identity of the PTV is lost and the sub-region can no longer be individually contained in an MLC hole.

[0184] If the overlap size between regions of different PTV groups exceeds a threshold whose minimum value is one pixel, a collapse strategy can reduce the number of regions. For spatially disjoint targets, the regions rarely overlap, so they do not provide any computational improvement. On the other hand, in the case of multiple enhancement volumes residing inside the base PTV, computational improvements may occur.

[0185] E. Trajectory Optimization Methods for Multi-Target Radiotherapy Treatment

[0186] Fig.16 A flow chart illustrating a method 1600 of trajectory optimization for radiation therapy treatment of multiple targets in accordance with some embodiments is shown.

[0187] At 1602, a patient model is provided. The patient model includes a plurality of planning target volumes (PTVs).

[0188] At 1604, a delivery coordinate space (DCS) is defined. The DCS includes a plurality of vertices. Each respective vertex defines a respective beam direction view (BEV) plane.

[0189] At 1606, for each respective PTV in the plurality of PTVs, for each respective BEV plane of the respective vertex, and for each respective pixel of the respective BEV plane, a dose of the respective PTV from a respective beamlet originating from the respective pixel and passing through the respective pixel of the respective BEV plane is evaluated. A respective BEV score for the respective pixel of the respective BEV plane of the respective PTV is also evaluated based at least in part on the dose of the respective PTV from the respective beamlet.

[0190] At 1608 , for each respective PTV in the plurality of PTVs, a PTV-specific threshold BEV score for the respective PTV is determined based at least in part on the BEV scores of the pixels of the BEV plane of the respective PTV.

[0191] At 1610, for each respective PTV in the plurality of PTVs and for each respective BEV plane of the respective vertex, one or more respective BEV regions of the respective PTV are determined by comparing each respective BEV score of the respective pixels of the respective PTV to a PTV-specific threshold BEV score. The BEV scores of the pixels within the one or more respective BEV regions are greater than or equal to the PTV-specific threshold BEV score.

[0192] At 1612, for each respective PTV in the plurality of PTVs, a respective BEV region connectivity manifold of the respective PTV is determined. The BEV region connectivity manifold represents connections between BEV regions of adjacent vertices.

[0193] At 1614, one or more optimal treatment trajectories for irradiating the plurality of PTVs are selected based on the BEV region connectivity manifolds of the plurality of PTVs.

[0194] You should understand that Fig.16 The specific steps shown provide a specific method for trajectory optimization for radiotherapy treatment according to some embodiments. Other sequences of steps may also be performed according to alternative embodiments. For example, alternative embodiments may perform the steps in different orders. Moreover, each step may include multiple sub-steps, which may be performed in various orders depending on the circumstances of the individual steps. Further, additional steps may be added and some steps may be removed depending on the specific application. Those of ordinary skill in the art will recognize many variations, modifications, and alternatives.

[0195] VI. Beam Angle Optimization Using Slices

[0196] The TORUS method uses BEV regions and BEV region connectivity manifolds as guidance to automatically generate heuristic optimal radiotherapy trajectories in order to efficiently deliver high-quality VMAT treatment plans. According to some embodiments, the TORUS method is modified to generate a heuristic optimal IMRT field. As discussed above, the TORUS method uses an optimization graph on top of the delivery coordinate space to generate an optimal treatment trajectory by using dual metric optimization. The nodes in the optimization graph represent individual control points, and the trajectory is defined as the path that minimizes the minimum distance metric, while the maximum distance metric acts as a goodness metric for selecting the optimal trajectory. For IMRT treatment plans, the treatment field can be represented in this framework as a set of nodes. Beam angle optimization may involve finding a set of k nodes with the best maximum distance.

[0197] A. Field Geometry Optimization

[0198] According to some embodiments, to find a set of k bundles, a BAO graph based on the TORUS graph concept may be constructed based on the following nodes and edges:

[0199] node

[0200] A node is a set of N tuples of the form (v, c, b) Where v is the vertex, c is the collimator angle index, and b is the binary mask of the included regions. Each tuple may identify a single field. A vertex v may correspond to a location in a discretized delivery coordinate space (DCS). For example, in a C-arm linear accelerator, a vertex v may correspond to the isocenter, the gantry angle θ, and the collimator angle θ. gantry and bed angle θ couch The collimator angle index c can be associated with a collimator angle θ in the set of possible discrete collimator angles. c The region mask b may correspond to a set of continuous target regions that the MLC blade may expose.

[0201] edge

[0202] Edges are connections between neighbor nodes that have vertex-vertex connectivity in the underlying delivery coordinate space and have MLC connectivity between corresponding MLC leaf sequences. In the case of beam angle optimization, the MLC connectivity constraint may not be important, but it can still be used to reduce the number of edges in the graph. Therefore, it may be computationally useful to maintain the MLC connectivity constraint.

[0203] In the TORUS method, minimum and maximum distances can be defined along the edges in the TORUS graph. In the case of a static field, when the treatment beam is on, there is no gantry motion. Therefore, only the vertices themselves need to be considered, and no minimum distance function is required. The score (maximum distance) can be the metric to be optimized. The score S can be defined on a subset of nodes as follows:

[0204]

[0205] The best set of k bundles can be defined as the subset that gives the best score S(B) in in Representation Subset The number of bundles in .

[0206] The score S(B) may be a complex non-local function of the entire subset B. Therefore, an exact solution may not be a trivial problem. However, there may be an efficient way to find an approximate solution faster than trying all NCk combinations. In some embodiments, beam angle optimization may include the following steps:

[0207] (1) Evaluate individual bundles and use the squares of these scores as sampling probabilities.

[0208] (2) Randomly sample k bundles using the sampling probability.

[0209] (3) Apply local gradient descent to the set B to find the local minimum.

[0210] (4) Repeat steps (1) to (3) until no improvement in the score is found in a predetermined number of consecutive trials (e.g., 20 consecutive trials).

[0211] This process may be referred to as a gradient descent method. During each iteration of the gradient descent process, all possible "edges" may be considered for improvement, where an "edge" is the movement of a subset B to a neighbor subset B'. The subsets of the bundle may be updated to move in the direction that will most improve the score. The process may be repeated until no more local improvements in the score are found. A neighbor subset of a subset B is defined as another subset B' that has exactly one different subset B in the bundle. i →B′ i , where the two bundles share an edge in the TORUS graph or have the same position (same vertex) in the delivery coordinate space.

[0212] B. Score Function

[0213] The score function S(B) can be similar to the maximum distance function in TORUS. The score S(B) of a bundle set B can include two parts: the individual bundle score and the overall PTV angular flux novelty. The existence of a global PTV angular flux metric makes the score function non-local, so that it is not just a function of the scores of individual bundles. The score S(B) can be written as the sum of the two parts,

[0214]

[0215] The local part It can be defined as follows:

[0216]

[0217]

[0218] w(B)=max(0.1,min(1.0,1.1-c 2 (B))), (20)

[0219] c(B) = MLC contention severity (21)

[0220] s(B) = integral score of the region (22)

[0221]

[0222] A avg can be defined as the average area of ​​the region. Therefore, the term can be considered as the average area score. In some embodiments, the MLC contention severity can be compared to A, which is defined above as the total non-area area exposed by the fitted MLC. penalty The cross section can be defined as the y extent of the open MLC leaf normalized to the average diameter of the area as a circle. rms (B) can be defined as the RMS x range of the open leaf. Flux section The PTV angular flux novelty can be defined as the angular flux state generated by the selected region, normalized by dividing by the average angular flux novelty of the individual regions. Splitting the score into local and non-local parts can allow some code optimizations to be done by pre-calculating the contribution of individual bundles to the local part.

[0223] C. Beam angle optimization including consideration of beam closing time

[0224] IMRT treatment is considered as external beam radiation therapy treatment, where the dose is delivered from k static beam positions r i =(v i , c i ), i=1,...,k delivery, where index v i and c i Corresponding to the apex and collimator angle index, respectively. In C-arm linear accelerator therapy systems, the apex can include the isocenter, gantry angle, and bed angle. In other types of external beam radiation therapy systems, the apex can include other treatment axis variables. From each r i The delivered fluence can be determined by a fluence optimization scheme that is based on the i To optimally distribute the dose between the two devices so as to meet the clinical optimization goals.

[0225] Each permissible combination can be Fluence optimization and leaf sorting are performed and a dosimetrically optimal IMRT treatment plan is found within a given delivery coordinate space by picking the dosimetrically optimal admissible combination. However, due to the large number of possible combinations, this process may be impractical. Moreover, it may be desirable for the patient and the patient's internal organs to remain stationary during the administration of radiation so that the delivered dose matches the planned dose. The longer the treatment takes, the greater the likelihood that the patient or the patient's internal organs will move during treatment, and therefore the greater the probability that the intended dose will not be delivered to the target volume.

[0226] The total treatment time can be calculated from the i To the beam i+1 The beam closing transition time is extended. i To the beam i+1 The total beam closing transition time can be expressed as Changes may include but are not limited to the following:

[0227] ··Move the machine axis from the vertex v i Move to the position v i+1 The time of the location

[0228] ··Move the collimator from c i Rotate to c i+1 Time

[0229] The time at which one or more images are acquired to guide treatment between beams i and i+1

[0230] Make the total bundle close time Minimum is likely to be an important part of the beam selection problem. In some embodiments, the BAO algorithm may include constraints on selecting beams that can be delivered in a time-sensitive radiation order. In some other embodiments, the radiation order of the beams may be determined as a post-processing step.

[0231] D. Beam Angle Optimization Method Using Slices

[0232] Fig.17 A flow chart illustrating a method 1700 of beam angle optimization using slices according to some embodiments is shown.

[0233] At 1702, a patient model is provided. The patient model includes one or more planning target volumes (PTVs), the one or more PTVs including a first PTV.

[0234] At 1704, a delivery coordinate space (DCS) is defined. The DCS includes a plurality of vertices. Each respective vertex defines a respective beam direction view (BEV) plane.

[0235] At 1706, for each respective BEV plane of the respective vertex and for each respective pixel of the respective BEV plane, a dose from a first PTV of a respective beamlet originating from the respective vertex and passing through the respective pixel of the respective BEV plane is evaluated. A respective BEV score for the respective pixel of the respective BEV plane is also evaluated based at least in part on the dose from the first PTV of the respective beamlet.

[0236] At 1708, a threshold BEV score for the first PTV is determined based on the BEV scores of the pixels of the BEV plane.

[0237] At 1710, for each respective BEV plane of the respective vertex, one or more BEV regions are determined by comparing the respective BEV score of each respective pixel to a threshold BEV score. The BEV scores of pixels within the one or more BEV regions are greater than or equal to the threshold BEV score.

[0238] At 1712, a BEV region connectivity manifold of the first PTV is determined. The BEV region connectivity manifold represents connections between BEV regions of adjacent vertices.

[0239] At 1714, a plurality of optimal field geometries for the IMRT radiation therapy treatment are selected based at least in part on the BEV region connectivity manifold of the first PTV.

[0240] You should understand that Fig.17 The specific steps shown provide a specific method for beam angle optimization according to some embodiments. Other sequences of steps may also be performed according to alternative embodiments. For example, alternative embodiments may perform the steps in a different order. Moreover, each step may include multiple sub-steps, which may be performed in various orders depending on the circumstances of the individual steps. Further, additional steps may be added and some steps may be removed depending on the specific application. Those of ordinary skill in the art will recognize many variations, modifications and alternatives.

[0241] E. Beam Angle Optimization Methods for Multi-Target Radiotherapy

[0242] Fig.18 A flow chart is shown illustrating a method 1800 of beam angle optimization for multi-target radiation therapy treatment in accordance with some embodiments.

[0243] At 1802, a patient model is provided. The patient model includes a plurality of planning target volumes (PTVs).

[0244] At 1804, a delivery coordinate space (DCS) is defined. The DCS includes a plurality of vertices. Each respective vertex defines a respective beam direction view (BEV) plane.

[0245] At 1806, for each respective PTV in the plurality of PTVs, for each respective BEV plane at the respective vertex, and for each respective pixel of the respective BEV plane, a dose of the respective PTV from a respective beamlet originating from the respective vertex and passing through the respective pixel of the respective BEV plane is evaluated. A respective BEV score for the respective pixel of the respective BEV plane of the respective PTV is also evaluated based at least in part on the dose of the respective PTV from the respective beamlet.

[0246] At 1808 , for each respective PTV in the plurality of PTVs, a PTV-specific threshold BEV score for the respective PTV is determined based at least in part on the BEV scores of the pixels of the BEV plane of the respective PTV.

[0247] At 1810, for each respective PTV in the plurality of PTVs and for each respective BEV plane of the respective vertex, one or more respective BEV regions of the respective PTV are determined by comparing each respective BEV score of the respective pixels of the respective PTV to a PTV-specific threshold BEV score. The BEV scores of the pixels within the one or more BEV regions are greater than or equal to the PTV-specific threshold BEV score.

[0248] At 1812, for each respective PTV in the plurality of PTVs, a respective BEV region connectivity manifold of the respective PTV is determined. The BEV region connectivity manifold represents connections between BEV regions of adjacent vertices.

[0249] At 1814, a plurality of optimal field geometries for IMRT radiation therapy treatment are selected based on the BEV region connectivity manifolds of the plurality of PTVs.

[0250] You should understand that Fig.18 The specific steps shown provide a specific method for beam angle optimization according to some embodiments. Other sequences of steps may also be performed according to alternative embodiments. For example, alternative embodiments may perform the steps in a different order. Moreover, each step may include multiple sub-steps, which may be performed in various orders depending on the circumstances of the individual steps. Further, additional steps may be added and some steps may be removed depending on the specific application. Those of ordinary skill in the art will recognize many variations, modifications and alternatives.

[0251] F. Beam Angle Optimization Method Using Slices with Coverage Enhancement

[0252] Fig.19A flow chart illustrating a method 1800 of beam angle optimization using slices with coverage enhancement in accordance with some embodiments is shown.

[0253] At 1902, a patient model is provided. The patient model includes one or more planning target volumes (PTVs), the one or more PTVs including a first PTV.

[0254] At 1904, a delivery coordinate space (DCS) is defined. The DCS may include a plurality of vertices. Each respective vertex defines a respective beam direction view (BEV) plane.

[0255] At 1906, for each respective BEV plane of the respective vertex, and for each respective pixel of the respective BEV plane, a dose from a first PTV of a respective beamlet originating from the respective vertex and passing through the respective pixel of the respective BEV plane is evaluated. A respective BEV score for the respective pixel of the respective BEV plane is also evaluated based at least in part on the dose from the first PTV of the respective beamlet.

[0256] At 1908, a global threshold BEV score for the first PTV is determined based on the BEV scores of the pixels of the BEV plane.

[0257] At 1910, for each respective BEV plane of the respective vertex, one or more initial BEV regions are determined by comparing the respective BEV score of each respective pixel to a global threshold BEV score. The BEV scores of the pixels within the one or more initial BEV regions are greater than or equal to the global threshold BEV score. Each respective initial BEV region may expose one or more voxels of the first PTV to a beamlet corresponding to the respective vertex.

[0258] At 1912 , for each respective voxel of the first PTV, a voxel-specific threshold BEV fraction is determined based on the initial BEV area of ​​the BEV plane such that the respective voxel is exposed by beamlets corresponding to a predetermined proportion of vertices out of a total number of vertices.

[0259] At 1914, for each respective BEV plane of the respective vertex and for each respective pixel of the respective BEV plane, a local threshold BEV score for the respective pixel is determined based on the voxel-specific threshold BEV score. One or more BEV regions of the respective BEV plane are determined by comparing the respective BEV score of each respective pixel to the local threshold BEV score. The BEV scores of the pixels within the one or more BEV regions are greater than or equal to the local threshold BEV score.

[0260] At 1916, a BEV region connectivity manifold of the first PTV is determined. The BEV region connectivity manifold represents connections between BEV regions of adjacent vertices.

[0261] At 1918, a plurality of optimal field geometries for the IMRT radiation therapy treatment are selected based at least in part on the BEV region connectivity manifold of the first PTV.

[0262] You should understand that Fig.19 The specific steps illustrated provide specific methods for beam angle optimization according to some embodiments. Other sequences of steps may also be performed according to alternative embodiments. For example, alternative embodiments may perform the steps in a different order. Moreover, each step may include multiple sub-steps, which may be performed in various orders depending on the circumstances of the individual steps. Further, additional steps may be added and some steps may be removed depending on the specific application. Those of ordinary skill in the art will recognize many variations, modifications, and alternatives.

[0263] VII. Computer System

[0264] Any computer systems mentioned herein may utilize any suitable number of subsystems. Fig. 20 An example of such a subsystem in computer system 2000 is shown. In some embodiments, the computer system includes a single computer device, where the subsystem can be a component of the computer device. In other embodiments, the computer system can include multiple computer devices, each of which is a subsystem with internal components.

[0265] Fig. 20 The subsystems shown are interconnected via a system bus 2075. Additional subsystems such as a printer 2074, a keyboard 2078, one or more storage devices 2079, a monitor 2076 coupled to a display adapter 2082, and the like are shown. Peripheral devices and input / output (I / O) devices coupled to the I / O controller 2071 can be connected to the computer system by any number of means known in the art, such as a serial port 2077. For example, the serial port 2077 or an external interface 2081 (e.g., Ethernet, Wi-Fi, etc.) can be used to connect the computer system 2000 to a wide area network such as the Internet, a mouse input device, or a scanner. The interconnection via the system bus 2075 allows the central processor 2073 to communicate with each subsystem and control the execution of instructions from the system memory 2072 or one or more storage devices 2079 (e.g., a fixed disk such as a hard drive or optical disk), as well as the exchange of information between the subsystems. The system memory 2072 and / or the one or more storage devices 2079 can be embodied as a computer-readable medium. Any data mentioned herein may be exported from one component to another component and may be exported to a user.

[0266] The computer system may include multiple identical components or subsystems, for example, connected together via external interface 2081 or via internal interfaces. In some embodiments, the computer system, subsystem or device may communicate via a network. In this case, one computer may be considered a client and another computer may be considered a server, where each computer may be part of the same computer system. The client and server may each include multiple systems, subsystems or components.

[0267] The external interface 2081 can be used to transmit one or more treatment plans to one or more radiation therapy devices, as described herein. For example, a treatment planning application can reside on a server computer, and a client computer can use the treatment planning application. The server computer can be part of a cloud computing platform that provides software as a service (SaaS). Once a treatment plan is determined, the client computer can specify which radiation device or radiation device accessible treatment plan database to use to transmit one or more files encapsulating the treatment plan. For example, an IP address can be specified.

[0268] It should be understood that any of the embodiments of the present invention can be implemented in a modular or integrated manner using hardware (e.g., an application specific integrated circuit or a field programmable gate array) and / or using computer software with a general programmable processor in the form of control logic. As used herein, the processor includes a multi-core processor on the same integrated chip or multiple processing units on a single circuit board or networked. Based on the disclosure and teachings provided herein, those of ordinary skill in the art will know and appreciate other ways and / or methods of implementing embodiments of the present invention using hardware and a combination of hardware and software.

[0269] Any of the software components or functions described in this application can be implemented as software code executed by a processor using any suitable computer language (such as, for example, Java, C++ or Perl) using, for example, conventional techniques or object-oriented techniques. The software code can be stored as a series of instructions or commands on a computer-readable medium for storage and / or transmission, and suitable media include random access memory (RAM), read-only memory (ROM), magnetic media such as a hard drive or floppy disk, or optical media such as a compact disk (CD) or DVD (digital versatile disk), flash memory, etc. The computer-readable medium can be any combination of such storage devices or transmission devices.

[0270] These programs can also be encoded and transmitted using carrier signals suitable for transmitting via wired networks, optical networks and / or wireless networks that meet multiple protocols, and these networks include the Internet. In this way, a computer-readable medium according to an embodiment of the present invention can be created using a data signal encoded by such a program. The computer-readable medium encoded by program code can be encapsulated with compatible devices or provided separately from other devices (for example, downloaded via the Internet). Any such computer-readable medium can reside on or in a single computer product (for example, a hard drive, a CD or a whole computer system), and may be present on or in different computer products in a system or network. A computer system can include a monitor, a printer or other suitable displays for providing any result in the result mentioned herein to a user.

[0271] Any method in the method described herein can be performed in whole or in part using a computer system including one or more processors, which can be configured to perform these steps. Therefore, an embodiment can be related to a computer system configured to perform the steps of any method in the method described herein, and the computer system may have different components for performing corresponding steps or corresponding step groups. Although presented in numbered steps, the steps of the method herein can be performed simultaneously or in different orders. Additionally, parts of these steps can be used together with parts of other steps from other methods. In addition, all or part of the steps can be optional. Additionally, any step of the steps of any method in the method can be performed using modules, circuits or other means for performing these steps.

[0272] The specific details of a particular embodiment may be combined in any suitable manner without departing from the spirit and scope of embodiments of the invention. However, other embodiments of the invention may be directed to specific embodiments relating to each individual aspect or specific combinations of these individual aspects.

[0273] The above description of exemplary embodiments of the present invention has been presented for purposes of illustration and description. It is not intended to be exhaustive or to limit the invention to the precise form described, and many modifications and variations are possible in light of the above teachings. The embodiments were chosen and described in order to best explain the principles of the invention and its practical application, thereby enabling others skilled in the art to best utilize the invention in various embodiments and with various modifications as are suitable for the particular use contemplated. .

[0274] Unless specifically indicated to the contrary, the term "a", "an" or "the" is intended to mean "one or more".

[0275] All patents, patent applications, publications, and descriptions mentioned herein are incorporated by reference in their entirety for all purposes. None are admitted to be prior art.

[0276] References

[0277] [1] Christopher Barry Locke and Karl Kenneth Bush. Trajectory optimization in radiotherapy using sectioning (TORUS). Medical Physics, 2017.

Claims

1. A trajectory optimization method for multi-target radiotherapy treatment, the method comprising: Providing a patient model, the patient model comprising a plurality of planning target volumes, i.e., a plurality of PTVs; defining a delivery coordinate space DCS, the DCS having a plurality of vertices, each respective vertex defining a respective beam direction view BEV plane, the delivery coordinate space being a set of all allowable coordinates parameterizing a configuration of a delivery device and being truncated to avoid collisions; For each respective PTV in the plurality of PTVs: For each corresponding BEV plane of the corresponding vertex, for each respective pixel of the respective BEV plane: evaluating a dose to the respective PTV from a respective beamlet originating from the respective vertex and passing through the respective pixel of the respective BEV plane; as well as evaluating a respective BEV score for each respective pixel of the respective BEV plane of the respective PTV based at least in part on the assessed dose of the respective PTV from the respective beamlet; determining a PTV-specific threshold BEV score for the corresponding PTV based at least in part on the BEV scores of the pixels of the BEV plane of the corresponding PTV; for each respective BEV plane of the respective vertices, determining one or more respective BEV regions of the respective PTV by comparing each respective BEV score of the respective pixels of the respective PTV to the PTV-specific threshold BEV score; and Determining a corresponding BEV region connectivity manifold of the corresponding PTV, wherein the corresponding BEV region connectivity manifold represents connections between BEV regions of adjacent vertices; as well as Based on the BEV region connectivity manifold of the plurality of PTVs, one or more optimal treatment trajectories are selected for irradiating the plurality of PTVs. 2 . The method of claim 1 , wherein the BEV scores of pixels within the one or more respective BEV regions are greater than or equal to the PTV-specific threshold BEV score.

3. The method of claim 1 , wherein selecting the one or more optimal treatment trajectories comprises: Perform statechart optimization.

4. The method of claim 3 , wherein the state graph optimization is performed on a graph defined by a plurality of nodes, each node being associated with a corresponding vertex in the DCS and with a state, the state being related to a set of PTV angular fluxes for a set of sampling points distributed in the plurality of PTVs, each PTV angular flux being related to the novelty of a direction vector of an incident beamlet passing through a closed surface centered at a corresponding sampling point of the set of sampling points. 5 . The method of claim 4 , wherein the set of sampling points is distributed among the plurality of PTVs according to weight factors associated with the plurality of PTVs. 6 . The method of claim 4 , wherein the set of sampling points are distributed among the plurality of PTVs such that a point density of each respective PTV is a decreasing function of a volume of the respective PTV.

7. The method of claim 3, wherein selecting the one or more optimal treatment trajectories comprises: selecting a plurality of candidate trajectories based on the BEV region connectivity manifolds of the plurality of PTVs, the plurality of candidate trajectories minimizing a minimum distance function; as well as The one or more best trajectories that maximize the maximum distance function are selected from the plurality of candidate trajectories.

8. The method according to claim 7, further comprising: For each PTV in the plurality of PTVs: For each corresponding BEV plane of the corresponding vertex: determining an area score for each respective pixel of the respective BEV plane based on the respective BEV scores of the respective pixels; and determining an average area score for each respective BEV region of the one or more BEV regions based on the area scores of the pixels of the respective BEV regions; Wherein the minimum distance function is inversely proportional to the average region score of each respective BEV region of each respective PTV, and the maximum distance function is directly proportional to the average region score of each respective BEV region of the respective PTV.

9. The method of claim 8, wherein the region score of each corresponding pixel has a normalized positive value if the corresponding pixel is within one of the one or more BEV regions, and has a negative value if the corresponding pixel is not within any of the one or more BEV regions.

10. The method of claim 1, wherein the patient model further includes organs at risk (OARs), the method further comprising: For each corresponding BEV plane of the corresponding vertex: For each corresponding pixel of the corresponding BEV plane: evaluating a dose to the OAR from a corresponding beamlet originating from the corresponding apex and passing through the corresponding pixel of the corresponding BEV plane; wherein evaluating the respective BEV score of the respective pixel of the respective BEV plane of the respective PTV comprises: evaluating a weighted combination of the dose to the respective PTV and the dose to the OAR from the respective beamlet; as well as The dose to the corresponding PTV is assigned a positive weight, and the dose to the OAR is assigned a negative weight.

11. A computer product comprising a non-transitory computer readable medium storing a plurality of instructions that, when executed, control a computer system to perform trajectory optimization for radiation therapy treatment of multiple targets, the instructions comprising: Providing a patient model, the patient model comprising a plurality of planning target volumes, i.e., a plurality of PTVs; defining a delivery coordinate space DCS, the DCS having a plurality of vertices, each respective vertex defining a respective beam direction view BEV plane; For each respective PTV in the plurality of PTVs: For each corresponding BEV plane of a corresponding vertex, for each corresponding pixel of said corresponding BEV plane: evaluating a dose to the respective PTV from a respective beamlet originating from the respective apex and passing through the respective pixel of the respective BEV plane; as well as evaluating a respective BEV score for each respective pixel of the respective BEV plane of the respective PTV based at least in part on the assessed dose of the respective PTV from the respective beamlet; determining a PTV-specific threshold BEV score for the corresponding PTV based at least in part on the BEV scores of the pixels of the BEV plane of the corresponding PTV; for each respective BEV plane of the respective vertices, determining one or more respective BEV regions of the respective PTV by comparing each respective BEV score of the respective pixels of the respective PTV to the PTV-specific threshold BEV score; and Determining a corresponding BEV region connectivity manifold of the corresponding PTV, wherein the corresponding BEV region connectivity manifold represents connections between BEV regions of adjacent vertices; as well as Based on the BEV region connectivity manifold of the plurality of PTVs, one or more optimal treatment trajectories are selected for irradiating the plurality of PTVs.

12. The computer product of claim 11, wherein pixels within the one or more respective BEV regions have BEV scores greater than or equal to the PTV-specific threshold BEV score.

13. The computer product of claim 11, wherein selecting the one or more optimal treatment trajectories comprises: Perform statechart optimization.

14. The computer product of claim 13, wherein the state graph optimization is performed on a graph defined by a plurality of nodes, each node being associated with a corresponding vertex in the DCS and with a state related to a set of PTV angular fluxes for a set of sampling points distributed in the plurality of PTVs, each PTV angular flux being related to the novelty of a direction vector of an incident beamlet passing through a closed surface centered at a corresponding sampling point of the set of sampling points.

15. The computer program product of claim 14, wherein the set of sampling points are distributed among the plurality of PTVs according to weight factors associated with the plurality of PTVs.

16. The computer product of claim 14, wherein the set of sampling points are distributed among the plurality of PTVs such that a point density of each respective PTV is a decreasing function of a volume of the respective PTV.

17. The computer product of claim 13, wherein selecting the one or more optimal treatment trajectories comprises: selecting a plurality of candidate trajectories based on the BEV region connectivity manifolds of the plurality of PTVs, the plurality of candidate trajectories minimizing a minimum distance function; as well as The one or more best trajectories that maximize the maximum distance function are selected from the plurality of candidate trajectories.

18. The computer program product of claim 17, wherein the instructions further comprise: For each PTV in the plurality of PTVs: For each corresponding BEV plane of the corresponding vertex: determining an area score for each respective pixel of the respective BEV plane based on the respective BEV scores of the respective pixels; and Based on the area scores of the pixels of the corresponding BEV regions, an average area score of each corresponding BEV region of the one or more BEV regions is determined; wherein the minimum distance function is inversely proportional to the average area score of each corresponding BEV region of each corresponding PTV, and the maximum distance function is directly proportional to the average area score of each corresponding BEV region of the corresponding PTV.

19. The computer product of claim 18, wherein the region score for each corresponding pixel has a normalized positive value if the corresponding pixel is within one of the one or more BEV regions, and has a negative value if the corresponding pixel is not within any of the one or more BEV regions.

20. The computer product of claim 11, wherein the patient model further includes organs at risk (OARs), the instructions further comprising: For each corresponding BEV plane of the corresponding vertex: for each respective pixel of the respective BEV plane: evaluating a dose to the OAR from a respective beamlet originating from the respective apex and passing through the respective pixel of the respective BEV plane; wherein evaluating the respective BEV score of the respective pixel of the respective BEV plane of the respective PTV comprises: evaluating a weighted combination of the dose to the respective PTV and the dose to the OAR from the respective beamlet; as well as The dose to the corresponding PTV is assigned a positive weight, and the dose to the OAR is assigned a negative weight.

Images