A radiotherapy plan generation method and system based on adaptive Monte Carlo calculation

Abstract

The application discloses a radiotherapy plan generation method and system based on adaptive Monte Carlo calculation. The method comprises the following steps: constructing an initial dose contribution matrix; performing analytical iterative optimization based on the initial dose contribution matrix; monitoring the convergence trend of the optimization state parameter in the process of analytical iterative optimization, and judging whether the analytical optimization enters a convergence stagnation state according to the convergence trend; when it is judged that the convergence stagnation state is entered, triggering Monte Carlo dose calculation, and constructing a Monte Carlo dose contribution matrix; correcting the initial dose contribution matrix according to the difference between the initial dose contribution matrix and the Monte Carlo dose contribution matrix, and generating a corrected dose contribution matrix; continuing to perform iterative optimization based on the corrected dose contribution matrix, and outputting a radiotherapy plan. The application triggers the Monte Carlo dose calculation and dynamically corrects the dose contribution matrix adaptively, thereby ensuring that the plan quality approximates the full Monte Carlo optimization level, greatly shortening the optimization time consumption, and improving the calculation efficiency and the plan quality.

Description

Technical Field

[0001] This invention relates to the field of radiotherapy technology, specifically to a method and system for generating radiotherapy plans based on adaptive Monte Carlo calculation. Background Technology

[0002] Intensity-modulated radiotherapy (IMRT) is currently the mainstream radiotherapy technique, and its planning quality highly depends on the accuracy of dose calculation during the optimization process. Existing commercial planning systems use analytical algorithms for dose calculation during the optimization phase, approximating human tissue as a layered homogeneous medium and performing only one-dimensional density correction along the radiation direction. Lateral scattering is approximated using a Gaussian distribution, resulting in extremely fast calculation speeds. However, in non-homogeneous tissue regions such as the lungs, nasopharynx, and bone-cavity junctions, the inherent physical approximations of analytical algorithms can introduce dose deviations of 5%–20%, particularly noticeable at small fields, built-up areas, and tissue interfaces.

[0003] The Monte Carlo (MC) algorithm, by randomly simulating all physical interactions between particles and the medium, can accurately handle tissue heterogeneity and lateral scattering imbalances, and is recognized as the "gold standard" for dose calculation. If MC is called for dose calculation in every iteration of intensity-modulated radiosampling (IMRT) optimization, the total optimization time will still increase dramatically from minutes to tens of hours, which cannot meet the real-time requirements of clinical practice.

[0004] In current clinical practice, "analytical optimization throughout + final MC validation" is a common approach. Using only analytical methods can prevent the optimizer from generating a true dose distribution. When MC validation shows underdose in the target area or overdose in organs at risk, the physicist must manually adjust and re-optimize and validate, which is inefficient and makes it difficult to guarantee plan quality. Another approach is to force a switch to MC calculation in the last few iterations of optimization. However, the switching timing is usually a pre-set fixed number of iterations, which cannot adapt to the optimization convergence process. Switching too early wastes computational resources, while switching too late introduces irreversible dose bias into the analytical calculation model. Furthermore, when jumping directly from the analytical engine to the MC engine, the significant differences in the response characteristics of the two dose calculation models often lead to oscillations or even divergence in the optimization objective function.

[0005] In summary, the existing technology has the following main drawbacks: 1. The pure analytical engine has insufficient optimization accuracy, which can easily produce clinically unacceptable dosage errors in complex anatomical sites; 2. Pure MC iterative optimization takes too long and lacks clinical real-time performance; 3. The strategy of switching MC at fixed times lacks adaptive awareness of the optimized state, resulting in poor efficiency and robustness. 4. The lack of a smooth transition mechanism when switching from the parsing engine to MC makes it difficult to guarantee stability.

[0006] In view of this, this invention patent is hereby proposed. Summary of the Invention

[0007] This invention aims to solve the above problems and proposes a radiotherapy planning generation method and system based on adaptive Monte Carlo calculation. Through the MC triggering mechanism and dynamic correction of the control point dose contribution matrix, the intensity-modulated radiotherapy (IMRT) optimization results can approach the accuracy of full MC optimization without significantly increasing the calculation time.

[0008] Specifically, the following technical solution was adopted: A method for generating radiotherapy plans based on adaptive Monte Carlo computation, comprising: Step S1: Acquire patient image data, establish geometric models of the target area and organs at risk, set prescription dose and dose constraints, use analytical dose algorithm to calculate the unit dose contribution of each sub-beam to the voxel, and construct the initial dose contribution matrix. Step S2: Based on the initial dose contribution matrix, with the weights of each sub-beam as optimization variables, an iterative optimization algorithm is used to perform analytical iterative optimization of the dose distribution. In each iteration, the cumulative dose is calculated based on the current weight vector and the initial dose contribution matrix, and the optimized state parameters are calculated based on the cumulative dose. Step S3: During the analytical iterative optimization process, monitor the convergence trend of the optimization state parameters, and determine whether the analytical optimization has entered a convergence stagnation state based on the convergence trend; when it is determined that it has entered a convergence stagnation state, trigger the Monte Carlo dose calculation; otherwise, return to step S2 to continue executing the analytical iterative optimization. Step S4: When Monte Carlo dose calculation is triggered, the Monte Carlo algorithm is used to calculate the unit dose contribution coefficient of each sub-beam to the voxel and construct the Monte Carlo dose contribution matrix; Step S5: Based on the difference between the initial dose contribution matrix and the Monte Carlo dose contribution matrix, the initial dose contribution matrix is ​​corrected to generate a corrected dose contribution matrix; Step S6: Based on the corrected dose contribution matrix, continue iterative optimization to obtain the optimal weight vector; Step S7: Output the radiotherapy plan based on the optimal weight vector.

[0009] Optionally, in the radiotherapy plan generation method based on adaptive Monte Carlo calculation of the present invention, the optimized state parameter in step S3 is the target optimization function value, and the convergence trend is the relative rate of change of the target optimization function value; The condition for determining whether a process has entered a convergence stagnation state is: the current iteration number t exceeds the preset minimum iteration number Nmin, and the relative rate of change of the objective optimization function value within p consecutive backtracking windows is [value missing]. All are below the preset threshold ; The relative rate of change The calculation formula is: ; In the formula, Let be the value of the objective function at the current iteration number t. Let m be the value of the target optimization function before the backtracking window m iterations, where m is the preset backtracking window step size; Optionally, the objective optimization function F(x) includes a weighted sum of the following constraint terms: Target area prescription dose constraints, target area dose uniformity constraints, maximum dose constraints for organs at risk, and dose penalty constraints for normal tissues; The general formula for the objective optimization function F(x) is: ; in: Target area prescription dose constraints: In the formula, For target area voxel collection, As the target area weighting factor, This represents the cumulative dose of voxel v. This is the prescription dosage; Target dose uniformity constraint: In the formula, It is a step function. This represents the maximum permissible dose to the target area. Maximum dose constraints for organs at risk: In the formula, To endanger organ somatic aggregates, As a factor that endangers organ weighting, Maximum permissible dose for organs at risk; Normal tissue dose penalty constraint: In the formula, It is a collection of normal tissue voxels. This is the weighting factor for normal tissues.

[0010] Optionally, in the radiotherapy plan generation method based on adaptive Monte Carlo calculation of the present invention, the optimized state parameters and the corresponding conditions for determining the convergence stagnation state in step S3 are any one of the following: The optimized state parameter is a clinical dosimetry indicator, and the convergence trend is the successive improvement of the clinical dosimetry indicator in multiple consecutive iterations. When the successive improvement is less than a preset clinical tolerance, it is determined to enter a convergence stagnation state. The clinical dosimetry indicator includes at least one of the target coverage indicator D95 or D98, the maximum dose to organs at risk Dmax, and the average dose to organs at risk Dmean; or, the optimized state parameter is the optimizer gradient norm. The convergence trend is the decay trend of the gradient norm. When the gradient norm drops below a preset threshold, it is determined that the convergence has stalled. Alternatively, the optimized state parameter is the sub-beam weight vector x, and the convergence trend is the successive change of the weight vector, where the change of the weight vector occurs in multiple consecutive iterations. When all values ​​are below the preset threshold, the system is determined to have entered a convergence stagnation state.

[0011] Optionally, in the radiotherapy planning generation method based on adaptive Monte Carlo calculation of the present invention, step S4, which involves calculating the unit dose contribution coefficient of each sub-beam to the voxel using the Monte Carlo algorithm, includes: Assign unit weights to each sub-beam j in sequence, reset the weights of the remaining sub-beams to zero, and perform Monte Carlo simulations to obtain the three-dimensional dose distribution contributed by each sub-beam individually. Extract the unit dose contribution coefficient of each sub-beam j to each voxel v. Construct the Monte Carlo dose contribution matrix Where V is the total number of voxels and N is the total number of sub-beams.

[0012] Optionally, in the radiotherapy planning generation method based on adaptive Monte Carlo calculation of the present invention, step S4, which involves calculating the unit dose contribution coefficient of each sub-beam to the voxel using the Monte Carlo algorithm, includes: Based on the current weight vector Perform a Monte Carlo simulation to obtain the synthetic dose distribution of each voxel v. ; Record particle trajectory information during Monte Carlo simulation; Based on the particle trajectory information and the current weight vector, the unit dose contribution coefficient of each sub-beam j to each voxel v is inferred through deconvolution or source model decomposition. Construct the Monte Carlo dose contribution matrix ; Wherein, the synthetic dose distribution Contribution coefficient of unit dose per sub-beam and weight The following relationship exists between them: ; By solving the aforementioned relationship, the unit dose contribution coefficient can be derived. .

[0013] Optionally, in the radiotherapy plan generation method based on adaptive Monte Carlo calculation of the present invention, step S5, which involves correcting the initial dose contribution matrix according to the difference between the initial dose contribution matrix and the Monte Carlo dose contribution matrix, includes: For each voxel v and each sub-beam j, the local dose correction factor is calculated as follows: : ; In the formula, The element values ​​in the initial dose contribution matrix, The element value at the corresponding position in the Monte Carlo dose contribution matrix, where δ is the preset low dose threshold; The local dose correction factor is spatially smoothed to obtain a smoothing correction factor. ; The corresponding elements of the initial dose contribution matrix are corrected point by point using the smoothing correction factor: ; Generate the corrected dose contribution matrix .

[0014] Optionally, in the radiotherapy plan generation method based on adaptive Monte Carlo calculation of the present invention, step S5 involves correcting the initial dose contribution matrix, including: A global modulation factor is introduced for each sub-beam j. The global modulation factor is solved by fitting the difference between the Monte Carlo dose contribution matrix and the initial dose contribution matrix using the least squares method. The corrected dose contribution matrix is ​​generated as follows: ; Alternatively, voxels can be divided into several groups according to tissue type, and the correction factor for each sub-beam j within each group can be the average of the correction factors of all voxels within that group. Correction dose contribution matrices can be generated according to tissue type; the tissue types include lung tissue, soft tissue, and bone tissue. Alternatively, a pre-trained machine learning model can be used to predict a correction factor field using the initial dose contribution matrix and local features of the patient image data as input. The initial dose contribution matrix can then be corrected based on the predicted correction factor field to generate a corrected dose contribution matrix.

[0015] Optionally, in the radiotherapy plan generation method based on adaptive Monte Carlo calculation of the present invention, step S6, which involves continuing iterative optimization based on the corrected dose contribution matrix until the final convergence condition is met, includes: The current weight vector when the Monte Carlo calculation is triggered As the initial value, the modified dose contribution matrix is ​​used. As the dose calculation engine, iterative optimization continues. In each iteration, the cumulative dose of each voxel v is calculated using the following formula. : ; In the formula, To correct the element values ​​of the dose contribution matrix, Let be the weight of the j-th sub-beam; Calculate the objective optimization function and its gradient based on the cumulative dose, and update the weight vector based on the gradient; When the objective function converges to the preset accuracy Or the number of iterations reaches the preset maximum number of iterations. When the iteration stops, the optimal weight vector is obtained; Optionally, after the final convergence condition is met, step S6 further includes: Independent Monte Carlo simulations were performed based on the obtained optimal weight vectors to obtain the validation dose distribution. The verified dose distribution is compared with the optimized dose distribution calculated based on the modified dose contribution matrix; When the deviation between the two exceeds the preset deviation threshold, the steps of calculating the Monte Carlo dose contribution matrix, matrix correction, and post-correction optimization iteration are repeatedly executed with the current optimal weight vector as the initial value, and the correction is performed twice or more in a loop until the deviation is lower than the preset deviation threshold.

[0016] Optionally, in the radiotherapy plan generation method based on adaptive Monte Carlo calculation of the present invention, the construction of the initial dose contribution matrix in step S1 includes: Initialize the beam sequence based on the prescribed dosage and target area geometry model; Each beam in the beam sequence is discretized into several sub-beam grids, and each sub-beam is an independent irradiation unit. The total number of sub-beams is N and the total number of voxels is V. The pencil beam analytical dosimetry algorithm is used to calculate the unit dose contribution coefficient of each sub-beam j to each voxel v. ; An initial dose contribution matrix is ​​constructed using the unit dose contribution coefficient of each sub-beam as elements. Each column corresponds to the unit dose contribution distribution of a sub-beam to all voxels, and each row corresponds to the unit dose contribution of a voxel to all sub-beams. The initial dose contribution matrix is ​​stored in the following sparse storage manner: a CSR format sparse index structure is established for the non-zero dose contribution relationship of each sub-beam to the voxel, and only the element values ​​exceeding the preset dose threshold are stored; The cumulative dose of voxel v under any sub-beam weight vector x Calculated using the following formula: .

[0017] Optionally, in the radiotherapy plan generation method based on adaptive Monte Carlo computation of the present invention, step S7, which involves outputting the radiotherapy plan according to the optimal weight vector, includes: Based on the optimal weight vector, the weight distribution of each sub-beam is mapped to the blade motion trajectory of the multi-leaf collimator using a blade sorting algorithm, generating an executable sub-field sequence. Using a Monte Carlo engine independent of step S4, the final dose verification is performed based on the optimal weight vector to obtain the deliverable dose distribution; The deliverable dose distribution is compared with the optimized dose distribution calculated based on the modified dose contribution matrix using dose-volume histogram indices and gamma analysis. When the gamma analysis pass rate meets the preset standard, a treatment plan report is generated; the treatment plan report includes a dose-volume histogram curve, a dose distribution map, a gamma analysis report, and the subfield sequence.

[0018] This invention also provides a radiotherapy planning generation system based on adaptive Monte Carlo computation, comprising: Acquisition Module: Used to acquire patient image data, establish geometric models of the target area and organs at risk, set prescription dose and dose constraints, use analytical dose algorithms to calculate the unit dose contribution of each sub-beam to the voxel, and construct the initial dose contribution matrix; The analytical optimization module is used to perform analytical iterative optimization of the dose distribution based on the initial dose contribution matrix, with the weights of each sub-beam as optimization variables, and using an iterative optimization algorithm. In each iteration, the cumulative dose is calculated based on the current weight vector and the initial dose contribution matrix, and the optimized state parameters are calculated based on the cumulative dose. Triggering Decision Module: Used to monitor the convergence trend of the optimization state parameters during the analytical iterative optimization process, and determine whether the analytical optimization has entered a convergence stagnation state based on the convergence trend; when it is determined that it has entered a convergence stagnation state, the Monte Carlo dose calculation is triggered; otherwise, the analytical optimization module is triggered to continue executing the analytical iterative optimization. Monte Carlo calculation module: When Monte Carlo dose calculation is triggered, it uses the Monte Carlo algorithm to calculate the unit dose contribution coefficient of each sub-beam to the voxel and constructs the Monte Carlo dose contribution matrix. Matrix correction module: used to correct the initial dose contribution matrix based on the difference between the initial dose contribution matrix and the Monte Carlo dose contribution matrix, and generate a corrected dose contribution matrix; Continue optimization module: used to continue iterative optimization based on the corrected dose contribution matrix to obtain the optimal weight vector; Output module: Used to output a radiotherapy plan based on the optimal weight vector.

[0019] Compared with existing technologies, this invention proposes a radiotherapy plan generation method and system based on adaptive Monte Carlo calculation, which has the following technical advantages: (1) Intelligent triggering and efficient use of computing resources: By judging whether the analytical optimization has entered a convergence stagnation state based on the convergence trend, the timing of MC intervention is adaptively determined, completely avoiding the resource waste or accuracy loss of fixed-step switching strategies. For simple cases, the analytical method converges quickly, and MC is triggered in advance to reduce unnecessary analytical iterations; for complex cases, the analytical method converges slowly, and MC is triggered accordingly with a delay to ensure that the analytical method searches sufficiently. MC is only called when the method engine truly needs accurate dose information, minimizing MC call overhead while ensuring the accuracy of the plan.

[0020] (2) Smooth transition of the dose model and strong optimization robustness: The dose contribution matrix of the sub-beams is dynamically corrected instead of directly replacing the entire dose engine. The correction factor asymptotically corrects the numerical deviation of the analytical matrix, preserving the sparse structure and gradient smoothing characteristics of the analytical matrix, thus achieving a soft transition between the analytical method and MC. The continuity of optimization is guaranteed, and the oscillation and divergence of the objective function caused by model jumps are completely avoided.

[0021] (3) Achieving both speed and accuracy: the corrected matrix It retains the structured and fast computation characteristics of analytical matrices (sparse matrix-vector multiplication, millisecond-level), while the dose values ​​are precisely calibrated by MC. Subsequent optimization iterations still run with the efficiency of analytical methods, but the dose calculation accuracy is close to that of MC. In typical cases, the total optimization time is only about 15-20 minutes longer than pure analytical optimization (for one MC matrix calculation), but the plan quality is comparable to that of full MC optimization (time reduction of more than 80%).

[0022] (4) High versatility and easy integration: This invention only requires adding two functional modules, MC trigger determination and matrix correction, to the existing intensity-modulated optimization framework, without changing the core iterative optimization algorithm logic of the optimizer. No major adjustments to the clinical workflow are required.

[0023] (5) Good scalability: The MC triggering and judgment mechanism does not depend on a specific MC engine or parsing engine and can be adapted to any parsing-MC combination. The matrix correction method can be further extended to multiple correction scenarios (cascaded correction). For extremely complex anatomical sites (such as multiple lung metastases with extensive calcification), accuracy can be further improved through 2-3 rounds of correction. The parallel architecture supports elastic scaling from single GPU to multi-GPU clusters, adapting to the hardware conditions of different levels of medical institutions. Attached Figure Description

[0024] Figure 1 A flowchart of a radiotherapy plan generation method based on adaptive Monte Carlo calculation according to an embodiment of the present invention; Figure 2 A flowchart illustrating the convergence of the objective function and MC triggering in a radiotherapy planning method based on adaptive Monte Carlo computation according to an embodiment of the present invention; Figure 3 This invention provides a flowchart of the dynamic correction process for the sub-beam-organ point dose contribution matrix in a radiotherapy planning generation method based on adaptive Monte Carlo calculation. Detailed Implementation

[0025] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments.

[0026] Therefore, the following detailed description of embodiments of the present invention is not intended to limit the scope of the claimed invention, but merely illustrates some embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0027] It should be noted that, unless otherwise specified, the embodiments and features and technical solutions in the present invention can be combined with each other.

[0028] It should be noted that similar labels and letters in the following figures indicate similar items. Therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures.

[0029] In the description of this invention, it should be noted that the terms "upper," "lower," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings, or the orientation or positional relationship commonly used when the product of this invention is in use, or the orientation or positional relationship commonly understood by those skilled in the art. These terms are only for the convenience of describing this invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation, and therefore should not be construed as a limitation of this invention. In addition, the terms "first," "second," etc., are only used to distinguish descriptions and should not be construed as indicating or implying relative importance.

[0030] like Figure 1As shown, the radiotherapy planning generation method based on adaptive Monte Carlo calculation provided by this invention includes the following steps: Steps 101 to 104 are the analytical optimization initialization and iteration stages; Step 105 is the convergence judgment node—if the relative rate of change of the objective function is lower than the threshold ε for p consecutive times and the number of iterations exceeds the minimum number of iterations Nmin, then proceed to step 106 to trigger MC calculation; otherwise, return to step 103 to continue analytical iteration; Steps 106 to 107 are the core innovative steps: MC calculation of the sub-beam deposition dose matrix. (Step 106), and calculate the correction factor and generate the correction matrix. (Step 107); Step 108: Continue optimization using the corrected matrix until final convergence; Step 109: Output the intensity modulation scheme.

[0031] The following is combined Figures 1 to 3 The above steps will be explained in detail.

[0032] 1. Data acquisition and construction of the initial dose contribution matrix.

[0033] Step 101: Acquire patient imaging data, establish geometric models of the target area and organs at risk, and set prescription doses and dose constraints.

[0034] Specifically, the patient's CT localization images (slice thickness 1-3mm) and target / organ-at-risk delineation RTStructs are first acquired. The prescribed dosage is then set by a physicist or the planning system. And dosage constraints for each organ at risk.

[0035] Step 102: Calculate the unit dose contribution of each sub-beam to the voxel using the analytical dose algorithm, and construct the initial dose contribution matrix.

[0036] like Figure 3 The initial dose contribution matrix is ​​shown above. It has a sparse structure. The specific construction method is as follows: (1) Initialize the beam sequence based on the prescribed dose and target area geometry model; (2) Discretize each beam in the beam sequence into several sub-beamlet grids, with each sub-beamlet serving as an independent irradiation unit. The irradiation field of each control point is discretized into a 2.5mm × 2.5mm sub-beamlet grid, and the total number of sub-beamlets is denoted as N, and the total number of voxels is denoted as V; (3) The pencil beam analytical dose algorithm (i.e., FSPB algorithm) is used for each beamlet. ) Calculate the unit MU deposition dose for each three-dimensional voxel v. ; (4) Construct an initial dose contribution matrix using the unit dose contribution coefficient of each sub-beam as elements. Each column corresponds to the unit dose contribution distribution of a sub-beam to all voxels, and each row corresponds to the unit dose contribution of a voxel to all sub-beams. (5) Use sparse representation technology for storage A CSR-formatted sparse index structure is established for the non-zero dose contribution relationship of each sub-beam to the voxel, storing only elements that exceed a preset dose threshold (e.g., 0.01% of the maximum dose).

[0037] In subsequent analytical optimization iterations, the cumulative dose of voxel v under any sub-beam weight vector x is... Calculated using the following formula: ; In the formula, Let be the MU weight of the j-th sub-beam.

[0038] 2. Analysis and iterative optimization.

[0039] Based on the initial dose contribution matrix, the weights of each sub-beam are used as optimization variables, and an iterative optimization algorithm is employed to perform analytical iterative optimization of the dose distribution. In each iteration, the cumulative dose is calculated based on the current weight vector and the initial dose contribution matrix, and the optimized state parameters are calculated based on the cumulative dose.

[0040] Specifically, given the sub-beam weight vector x, the absorbed dose of any voxel v is quickly calculated using the dose contribution matrix. An objective optimization function F(x) is constructed, and the objective function is optimized using a weighted sum of the following constraint terms: (1) Target area prescription dose constraints: In the formula, For target area voxel collection, As the target area weighting factor, This represents the cumulative dose of voxel v. This is the prescription dosage.

[0041] (2) Target area dose uniformity constraint: In the formula, It is a step function. This is the maximum permissible dose for the target area.

[0042] (3) Maximum dose constraints for organs at risk: In the formula, To endanger organ somatic aggregates, As a factor that endangers organ weighting, Maximum permissible dose for organs at risk; (4) Normal tissue dose penalty constraint: In the formula, It is a collection of normal tissue voxels. This is the weighting factor for normal tissues.

[0043] The general form of the objective function is: .

[0044] The L-BFGS (Limited-memory Broyden-Fletcher-Goldfarb-Shanno) algorithm is used for optimization. In each iteration, the optimizer calculates the gradient of the objective function with respect to the sub-beam weights and updates the weight vector x according to the gradient direction. L-BFGS uses historical gradient information to approximate the inverse of the Hessian matrix, achieving superlinear convergence speed while maintaining low memory usage.

[0045] 3. Adaptive trigger determination.

[0046] During the analytical iterative optimization process, the convergence trend of the optimization state parameters is monitored, and it is determined whether the analytical optimization has entered a convergence stagnation state based on the convergence trend. When it is determined that the convergence stagnation state has been entered, the Monte Carlo dose calculation is triggered; otherwise, the analytical iterative optimization is returned to continue.

[0047] like Figure 2 As shown, the upper part is the analytical optimization stage; the middle part is the convergence detection region—the relative rate of change of the objective function in three consecutive backtracking windows (each with an interval of m=5 iterations). All below the threshold At that point, the system determines that the analytical method has converged to its accuracy limit, triggering MC optimization.

[0048] In this embodiment, the optimization state parameter is the target optimization function value, and the convergence trend is the relative rate of change of the target optimization function value. Every m iterations (m=5 recommended), the current target function value is recorded. , where t is the current iteration number. Calculate the relative rate of change of the objective function: ; In the formula, Let be the value of the objective function at the current iteration number t. This represents the target optimization function value before the backtracking window reaches m iterations.

[0049] The MC trigger condition must meet the following two conditions simultaneously: Condition 1: The current iteration number t exceeds the preset minimum iteration number Nmin (N_min=20 is recommended). This condition avoids accidentally triggering MC during the rapid decrease of the objective function in the early stage of optimization—in the early stage of optimization, the objective function may decrease significantly but the rate of change fluctuates wildly. At this time, there is still sufficient optimization space in the initial iteration, and there is no need for MC intervention.

[0050] Condition 2: Satisfies this condition for p consecutive times (p=3 recommended). ,in For convergence determination threshold (e.g.) This condition means that the initial iterative optimizer has failed to significantly improve the objective function for several consecutive backtracking windows and has entered a convergence plateau. The dose modeling accuracy of the analytical engine has become a bottleneck limiting further optimization.

[0051] False trigger protection mechanism: To prevent random fluctuations in the objective function from triggering MC (Match-Cut), a continuous multiple-judgment strategy is adopted (rather than a single threshold judgment). Simultaneously, a hard lower limit of the minimum iteration number Nmin is introduced to ensure that the analytical computation model performs sufficient initial search. Furthermore, it can assist in monitoring the independent changing trends of key clinical indicators (such as the successive improvement in PTV D95) as a supplementary reference for trigger decision-making.

[0052] In other embodiments of the present invention, the optimized state parameters and the conditions for determining the convergence stagnation state can also be implemented in the following ways: (1) The optimized state parameter is a clinical dosimetric index, and the convergence trend is the successive improvement of the clinical dosimetric index in multiple consecutive iterations. When the successive improvement is less than the preset clinical tolerance, it is determined to enter the convergence stagnation state. The clinical dosimetric index includes at least one of the target coverage index D95 or D98, the maximum dose to organs at risk Dmax, and the average dose to organs at risk Dmean. (2) The optimized state parameter is the optimizer gradient norm. The convergence trend is the decay trend of the gradient norm. When the gradient norm drops below a preset threshold, it is determined that the convergence has stalled. (3) The optimized state parameter is the sub-beam weight vector x, and the convergence trend is the successive change of the weight vector. When the change of the weight vector in consecutive iterations... When all values ​​are below the preset threshold, the system is determined to have entered a convergence stagnation state.

[0053] When the judgment condition is met, the system automatically triggers the MC dose calculation module; otherwise, it returns and continues the initial matrix optimization iteration.

[0054] 4. Calculation of MC dose contribution matrix.

[0055] When Monte Carlo dose calculation is triggered, the Monte Carlo algorithm is used to calculate the unit dose contribution coefficient of each sub-beam grid to each voxel based on the sub-beam weight distribution corresponding to the current optimal weight vector, and the Monte Carlo dose contribution matrix is ​​constructed.

[0056] like Figure 3 As shown, the MC module obtains the precise deposition dose coefficient of each sub-beam to the voxel by sequentially activating the sub-beams or by deconvolution of the source model. .

[0057] In this embodiment, one of the following two strategies is used to obtain the precise deposition dose coefficient of each sub-beam to the voxel by MC. : Strategy A (One-by-One Activation): Assign unit weights to each sub-beam j sequentially (reset the weights of other sub-beams to zero), and perform McLeod simulations to obtain the 3D dose distribution contributed individually by that sub-beam. This method offers the highest accuracy but also the highest computational cost, making it suitable for scenarios with a small number of sub-beams (N<1000). Extract the unit dose contribution coefficient of each sub-beam j to each voxel v. Construct the Monte Carlo dose contribution matrix .

[0058] Strategy B (Source Model Deconvolution): Based on the current weight vector Perform a Monte Carlo simulation to obtain the synthetic dose distribution of each voxel v. During the Monte Carlo simulation, particle trajectory information (such as the initial interaction position and energy deposition position) is recorded. Based on the particle trajectory information and the current weight vector, the unit dose contribution coefficient of each sub-beam j to each voxel is inferred through deconvolution or source model decomposition. Construct the Monte Carlo dose contribution matrix Among them, the synthetic dose distribution Contribution coefficient of unit dose per sub-beam and weight The following relationship exists between them: ; By solving the aforementioned relationship, the unit dose contribution coefficient can be derived. This method obtains the dose contribution of all sub-beams simultaneously in a single MC simulation, which is highly efficient and suitable for typical clinical scenarios with a large number of sub-beams (N≥1000).

[0059] The MC particle transport simulation uses a GPU-accelerated MC dose calculation engine. Recommended simulation parameters: 1 × 10⁻⁶ particle history. 7 (It can be decomposed into multiple independent subtasks and executed in parallel on a GPU cluster), with statistical uncertainty controlled below 1%.

[0060] 5. Dynamic correction of dose contribution matrix.

[0061] like Figure 3 As shown, the initial dose contribution matrix is ​​corrected based on the difference between the initial dose contribution matrix and the Monte Carlo dose contribution matrix to generate a corrected dose contribution matrix.

[0062] This invention does not use directly replace (To avoid abrupt model changes), a gradual matrix dynamic correction method is adopted, which includes the following sub-steps: (1) Calculation of local dose correction factor: For each voxel v and each sub-beam j, the local dose correction factor is calculated as follows: ; In the formula, The element values ​​in the initial dose contribution matrix, The element value at the corresponding position in the Monte Carlo dose contribution matrix is ​​δ, which is a preset low dose threshold (e.g., 0.1% of the maximum deposition dose) to avoid noise amplification in the low dose region.

[0063] (2) Spatial smoothing and regularization of the correction factor field: For the initial correction factor field Spatial smoothing filtering is performed to eliminate high-frequency spike noise caused by MC statistical fluctuations. One of the following smoothing schemes can be used: 3D Gaussian filtering: Perform 3D Gaussian kernel convolution on the spatial distribution of the V repair positive factor of each sub-beam j, with a kernel width of... It is recommended to set it to 1-2 times the voxel size (approximately 2-5 mm); Three-dimensional median filtering: Median filtering of the correction factor field (kernel size recommended as 3×3×3) can effectively eliminate isolated noise points while better preserving the correction factor gradient at the tissue interface; Adaptive bilateral filtering: It takes into account both spatial proximity and numerical similarity, and maintains the discontinuity of the correction factor at the tissue boundary while smoothing noise.

[0064] The correction factor field is obtained after smoothing. .

[0065] (3) Generate the corrected dose contribution matrix: The initial matrix is ​​corrected element by element using a smoothed correction factor: ; This generates the corrected dose contribution matrix. .like Figure 3 As shown below, The sparse structure of the analytical method is preserved (the position of the zero element remains unchanged), but the values ​​are precisely calibrated by MC.

[0066] In other embodiments of the present invention, the modification may also be carried out in the following manner: (1) Introduce a global modulation factor for each sub-beam j. The global modulation factor is solved by fitting the difference between the Monte Carlo dose contribution matrix and the initial dose contribution matrix using the least squares method. The corrected dose contribution matrix is ​​generated as follows: ; (2) The voxels are divided into several groups according to tissue type. The correction factor for each sub-beam j in each group is the average value of the correction factors of all voxels in that group. The correction dose contribution matrix is ​​generated according to tissue type. The tissue types include lung tissue, soft tissue, and bone tissue. (3) Using a pre-trained machine learning model, with the initial dose contribution matrix and local features of the patient image data as input, a correction factor field is predicted, and the initial dose contribution matrix is ​​corrected based on the predicted correction factor field to generate a corrected dose contribution matrix.

[0067] Will Send it back to the optimizer and replace the original one. This serves as the basis for dose calculations in subsequent optimizations. In subsequent optimization iterations, voxel doses are calculated using the corrected matrix: ; because Maintained with With the exact same sparse structure, the dose calculation efficiency remains unchanged in a single iteration; however, since each element has been precisely calibrated by MC, the calculated dose distribution accuracy approaches the MC level.

[0068] 6. Optimize and iterate after correction.

[0069] The current weight vector when the Monte Carlo calculation is triggered As the initial value, the modified dose contribution matrix is ​​used. As the dose calculation engine, it continues to perform iterative optimization.

[0070] like Figure 2 As shown in the second half, in the optimization stage after matrix correction, the objective function continues to decrease until it fully converges, and no oscillations caused by model switching occur throughout the process.

[0071] In each iteration, the cumulative dose of each voxel v is calculated using the following formula: ; The objective optimization function and its gradient are calculated based on the cumulative dose, and the weight vector is updated based on the gradient.

[0072] During the post-correction optimization phase, changes in the objective function value continue to be monitored. Because... The dosage accuracy has been greatly improved, and the objective function usually has a new second convergence process of decline, which further improves target coverage and OAR protection.

[0073] The revised optimization allows setting a final convergence criterion: when the objective optimization function converges to a preset accuracy. (in For example, take (or the number of iterations reaches the preset maximum number of iterations) (like When the condition is met, stop the iteration and obtain the optimal weight vector.

[0074] Optionally, when the corrected optimization is close to convergence, a second MC verification can be performed—an independent MC simulation based on the final sub-beam weights, comparing the synthetic dose with... Calculate the dosage and compare the results. If a significant deviation still exists (e.g., DVH index deviation >1%), repeat the MC calculation and matrix correction, performing secondary or multiple cyclic corrections to further improve accuracy. Clinical trial data indicate that a single correction is sufficient in the vast majority of cases.

[0075] 7. Output the treatment plan.

[0076] The radiotherapy plan is output based on the optimal weight vector, specifically including: (1) Based on the optimal weight vector, the weight distribution of each sub-beam is mapped to the blade motion trajectory of the multi-leaf collimator by the blade sorting algorithm, and an executable sub-field sequence is generated. (2) Using an independent Monte Carlo engine, the final dose verification is performed based on the optimal weight vector to obtain the deliverable dose distribution; (3) Compare the dose distribution with the optimized dose distribution calculated based on the modified dose contribution matrix using dose volume histogram index and gamma analysis (recommended standard: 3% / 3mm, 10% threshold). (4) When the gamma analysis pass rate meets the preset standard, a treatment plan report is generated; the treatment plan report includes a dose-volume histogram curve, a dose distribution map, a gamma analysis report and the subfield sequence, and is transmitted to the treatment record and verification system for clinical use.

[0077] The technical effects of this invention will be illustrated below using a case of central lung cancer in the right lung treated with VMAT (Vibration and Acute Thromboplasty) as an example.

[0078] The prescribed dose is 60 Gy / 30 fractions. The target volume (PTV) is approximately 185 cc, and organs at risk include the ipsilateral lung, contralateral lung, heart, esophagus, spinal cord, and proximal bronchial tree. Three treatment plans were designed for comparison: Plan A (Standard Clinical Mode): Full analysis and optimization, followed by delivery to an independent MC engine for validation dose calculation after optimization, representing the current standard clinical workflow; Plan B (ideal reference benchmark): The dose contribution matrix and synthetic dose are calculated using the MC engine in each iteration, and full MC optimization is performed until convergence, representing the theoretically achievable optimal plan quality; Plan C (the method of this invention): The initial matrix is ​​optimized to adaptively trigger MC determination, the MC dose contribution matrix is ​​calculated, the matrix is ​​dynamically corrected, and after correction, optimization continues until final convergence.

[0079] Key parameters are set as follows: , , , Spatial smoothing employs three-dimensional Gaussian filtering ( ), MC simulation particle number .

[0080] Table 1. Comparison of dosimetric parameters and time consumption for the three optimization methods:

[0081] The results showed that the target coverage and organ-at-risk protection levels of the method of this invention (Plan C) were highly similar to those of the full MC optimization (Plan B), with a PTV D98% of 59.1 Gy (59.2 Gy in full MC optimization), while the pure analytical optimization was only 56.3 Gy; the ipsilateral lung V20 was 19.1% (18.9% in full MC optimization), and the average cardiac dose was 11.4 Gy (11.2 Gy in full MC optimization). Regarding the total optimization time, the method of this invention took 22 minutes, approximately 84% shorter than the 135 minutes of full MC optimization, and only about 18 minutes longer than the 4 minutes of pure analytical optimization.

[0082] Table 2. Comparison of Final Dose and Deliver Dose Consistency (Plan C):

[0083] Gamma analysis was performed using a 3% / 3 mm standard and a 10% dose threshold. Deviations for all DVH indicators were less than 0.5%, and the global gamma pass rate reached 99.1%, demonstrating that the optimization results based on beamlet deposition dose dynamic correction are highly consistent with the actual deliverable dose, meeting the requirements for clinical dose validation.

[0084] The above description is merely a preferred embodiment of the present invention and does not limit the scope of the patent. In other alternative embodiments of the present invention: (1) The analytical dose algorithm is not limited to the FSPB algorithm, but can also use convolution algorithm or other analytical dose algorithms; (2) The iterative optimization algorithm is not limited to the L-BFGS algorithm, but can also use gradient descent, conjugate gradient or other iterative optimization algorithms; (3) The spatial smoothing process is not limited to three-dimensional Gaussian filtering, three-dimensional median filtering or adaptive bilateral filtering, and other spatial filtering methods may also be used; (4) The machine learning model is not limited to 3D U-Net or random forest, but may also use other deep learning or machine learning models; (5) The MC calculation can be implemented on a single GPU, multiple GPUs or CPU clusters, and the method of the present invention can be implemented under different hardware configurations.

[0085] This invention also provides a radiotherapy planning generation system based on adaptive Monte Carlo calculation, comprising: Acquisition Module: Used to acquire patient image data, establish geometric models of the target area and organs at risk, set prescription dose and dose constraints, use analytical dose algorithms to calculate the unit dose contribution of each sub-beam to the voxel, and construct the initial dose contribution matrix; The analytical optimization module is used to perform analytical iterative optimization of the dose distribution based on the initial dose contribution matrix, with the weights of each sub-beam as optimization variables, and using an iterative optimization algorithm. In each iteration, the cumulative dose is calculated based on the current weight vector and the initial dose contribution matrix, and the optimized state parameters are calculated based on the cumulative dose. Triggering Decision Module: Used to monitor the convergence trend of the optimization state parameters during the analytical iterative optimization process, and determine whether the analytical optimization has entered a convergence stagnation state based on the convergence trend; when it is determined that it has entered a convergence stagnation state, the Monte Carlo dose calculation is triggered; otherwise, the analytical optimization module is triggered to continue executing the analytical iterative optimization. Monte Carlo calculation module: When Monte Carlo dose calculation is triggered, it uses the Monte Carlo algorithm to calculate the unit dose contribution coefficient of each sub-beam to the voxel and constructs the Monte Carlo dose contribution matrix. Matrix correction module: used to correct the initial dose contribution matrix based on the difference between the initial dose contribution matrix and the Monte Carlo dose contribution matrix, and generate a corrected dose contribution matrix; Continue optimization module: used to continue iterative optimization based on the corrected dose contribution matrix to obtain the optimal weight vector; Output module: Used to output a radiotherapy plan based on the optimal weight vector.

[0086] The aforementioned modules can be integrated into a single processor, exist as separate physical units, or be integrated into a single unit. The integrated modules can be implemented in hardware or as software functional modules. When the integrated modules are implemented as software functional modules and sold or used as independent products, they can be stored in a computer-readable storage medium.

[0087] This embodiment also provides a computer-readable storage medium storing a computer-executable program, which, when executed, implements the radiotherapy plan generation method based on adaptive Monte Carlo calculation.

[0088] The computer-readable storage medium described in this embodiment may include data signals propagated in baseband or as part of a carrier wave, carrying readable program code. Such propagated data signals may take various forms, including but not limited to electromagnetic signals, optical signals, or any suitable combination thereof. The computer-readable storage medium may also be any readable medium other than a readable storage medium, capable of transmitting, propagating, or transmitting programs for use by or in connection with an instruction execution system, apparatus, or device. The program code contained on the computer-readable storage medium may be transmitted using any suitable medium, including but not limited to wireless, wired, optical fiber, RF, etc., or any suitable combination thereof.

[0089] This embodiment also provides an electronic device, including a processor and a memory, wherein the memory is used to store a computer-executable program, and when the computer program is executed by the processor, the processor executes the radiotherapy plan generation method based on adaptive Monte Carlo calculation.

[0090] The electronic device is manifested in the form of a general-purpose computing device. It may contain one or more processors that work collaboratively. This invention also does not preclude distributed processing, meaning that processors may be distributed across different physical devices. The electronic device of this invention is not limited to a single entity, but may also be the sum of multiple physical devices.

[0091] The memory stores a computer-executable program, typically machine-readable code. The computer-readable program can be executed by the processor to enable the electronic device to perform the method of the present invention, or at least some steps of the method.

[0092] The memory includes volatile memory, such as random access memory (RAM) and / or cache memory, and may also be non-volatile memory, such as read-only memory (ROM).

[0093] It should be understood that the electronic device of the present invention may also include elements or components not shown in the examples above. For example, some electronic devices also include display units such as a display screen, and some electronic devices also include human-computer interaction elements such as buttons and keyboards. Any electronic device capable of executing a computer-readable program in its memory to implement the method of the present invention or at least some steps of the method can be considered as an electronic device covered by the present invention.

[0094] From the above description of the embodiments, those skilled in the art will readily understand that the present invention can be implemented by hardware capable of executing specific computer programs, such as the system of the present invention, and the electronic processing unit, server, client, mobile phone, control unit, processor, etc. included in the system. The present invention can also be implemented by computer software executing the methods of the present invention, for example, by control software executed by a microprocessor, electronic control unit, client, server, etc. However, it should be noted that the computer software executing the methods of the present invention is not limited to execution in one or a specific set of hardware entities; it can also be implemented in a distributed manner by unspecified hardware. For computer software, the software product can be stored in a computer-readable storage medium (such as a CD-ROM, USB flash drive, portable hard drive, etc.) or distributed across a network, as long as it enables electronic devices to execute the methods according to the present invention.

[0095] The above embodiments are only used to illustrate the present invention and are not intended to limit the technical solutions described herein. Although the present invention has been described in detail with reference to the above embodiments, the present invention is not limited to the specific embodiments described above. Therefore, any modifications or equivalent substitutions to the present invention, as well as all technical solutions and improvements that do not depart from the spirit and scope of the invention, are covered within the scope of the claims of the present invention.

Claims

1. A method for generating radiotherapy plans based on adaptive Monte Carlo calculation, characterized in that, include: Step S1: Acquire patient image data, establish geometric models of the target area and organs at risk, set prescription dose and dose constraints, use analytical dose algorithm to calculate the unit dose contribution of each sub-beam to the voxel, and construct the initial dose contribution matrix. Step S2: Based on the initial dose contribution matrix, with the weights of each sub-beam as optimization variables, an iterative optimization algorithm is used to perform analytical iterative optimization of the dose distribution. In each iteration, the cumulative dose is calculated based on the current weight vector and the initial dose contribution matrix, and the optimized state parameters are calculated based on the cumulative dose. Step S3: During the analytical iterative optimization process, monitor the convergence trend of the optimization state parameters, and determine whether the analytical optimization has entered a convergence stagnation state based on the convergence trend; when it is determined that it has entered a convergence stagnation state, trigger the Monte Carlo dose calculation; otherwise, return to step S2 to continue executing the analytical iterative optimization. Step S4: When Monte Carlo dose calculation is triggered, the Monte Carlo algorithm is used to calculate the unit dose contribution coefficient of each sub-beam to the voxel and construct the Monte Carlo dose contribution matrix; Step S5: Based on the difference between the initial dose contribution matrix and the Monte Carlo dose contribution matrix, the initial dose contribution matrix is ​​corrected to generate a corrected dose contribution matrix; Step S6: Based on the corrected dose contribution matrix, continue iterative optimization to obtain the optimal weight vector; Step S7: Output the radiotherapy plan based on the optimal weight vector.

2. The radiotherapy plan generation method based on adaptive Monte Carlo calculation according to claim 1, characterized in that, The optimization state parameter mentioned in step S3 is the target optimization function value, and the convergence trend is the relative rate of change of the target optimization function value; The condition for determining whether a process has entered a convergence stagnation state is: the current iteration number t exceeds the preset minimum iteration number Nmin, and the relative rate of change of the objective optimization function value within p consecutive backtracking windows is [value missing]. All are below the preset threshold ; The relative rate of change The calculation formula is: ； In the formula, Let be the value of the objective function at the current iteration number t. Let m be the value of the target optimization function before the backtracking window m iterations, where m is the preset backtracking window step size; Optionally, the objective optimization function F(x) includes a weighted sum of the following constraint terms: Target area prescription dose constraints, target area dose uniformity constraints, maximum dose constraints for organs at risk, and dose penalty constraints for normal tissues; The general formula for the objective optimization function F(x) is: ； in: Target area prescription dose constraints: In the formula, For target area voxel collection, As the target area weighting factor, This represents the cumulative dose of voxel v. This is the prescription dosage; Target dose uniformity constraint: In the formula, It is a step function. This represents the maximum permissible dose to the target area. Maximum dose constraints for organs at risk: In the formula, To endanger organ somatic aggregates, As a factor that endangers organ weighting, Maximum permissible dose for organs at risk; Normal tissue dose penalty constraint: In the formula, It is a collection of normal tissue voxels. This is the weighting factor for normal tissues.

3. The radiotherapy plan generation method based on adaptive Monte Carlo calculation according to claim 1, characterized in that, The optimized state parameters and the corresponding conditions for entering the convergence stagnation state mentioned in step S3 are any one of the following: The optimized state parameter is a clinical dosimetry indicator, and the convergence trend is the successive improvement of the clinical dosimetry indicator in multiple consecutive iterations. When the successive improvement is less than a preset clinical tolerance, it is determined to enter a convergence stagnation state. The clinical dosimetry indicator includes at least one of the target coverage indicator D95 or D98, the maximum dose to organs at risk Dmax, and the average dose to organs at risk Dmean; or, the optimized state parameter is the optimizer gradient norm. The convergence trend is the decay trend of the gradient norm. When the gradient norm drops below a preset threshold, it is determined that the convergence has stalled. Alternatively, the optimized state parameter is the sub-beam weight vector x, and the convergence trend is the successive change of the weight vector, where the change of the weight vector occurs in multiple consecutive iterations. When all values ​​are below the preset threshold, the system is determined to have entered a convergence stagnation state.

4. A method for generating radiotherapy plans based on adaptive Monte Carlo calculation according to claim 2 or 3, characterized in that, Step S4, which involves using the Monte Carlo algorithm to calculate the unit dose contribution coefficient of each sub-beam to the voxel, includes: Assign unit weights to each sub-beam j in sequence, reset the weights of the remaining sub-beams to zero, and perform Monte Carlo simulations to obtain the three-dimensional dose distribution contributed by each sub-beam individually. Extract the unit dose contribution coefficient of each sub-beam j to each voxel v. Construct the Monte Carlo dose contribution matrix Where V is the total number of voxels and N is the total number of sub-beams.

5. A method for generating radiotherapy plans based on adaptive Monte Carlo calculation according to claim 2 or 3, characterized in that, Step S4, which involves using the Monte Carlo algorithm to calculate the unit dose contribution coefficient of each sub-beam to the voxel, includes: Based on the current weight vector Perform a Monte Carlo simulation to obtain the synthetic dose distribution of each voxel v. ; Record particle trajectory information during Monte Carlo simulation; Based on the particle trajectory information and the current weight vector, the unit dose contribution coefficient of each sub-beam j to each voxel v is inferred through deconvolution or source model decomposition. Construct the Monte Carlo dose contribution matrix ; Wherein, the synthetic dose distribution Contribution coefficient of unit dose per sub-beam and weight The following relationship exists between them: ； By solving the aforementioned relationship, the unit dose contribution coefficient can be derived. .

6. A method for generating radiotherapy plans based on adaptive Monte Carlo calculation according to claim 2 or 3, characterized in that, Step S5, which involves correcting the initial dose contribution matrix based on the difference between the initial dose contribution matrix and the Monte Carlo dose contribution matrix, includes: For each voxel v and each sub-beam j, the local dose correction factor is calculated as follows: : ； In the formula, The element values ​​in the initial dose contribution matrix, The element value at the corresponding position in the Monte Carlo dose contribution matrix, where δ is the preset low dose threshold; The local dose correction factor is spatially smoothed to obtain a smoothing correction factor. ; The corresponding elements of the initial dose contribution matrix are corrected point by point using the smoothing correction factor: ； Generate the corrected dose contribution matrix .

7. A method for generating radiotherapy plans based on adaptive Monte Carlo calculation according to claim 2 or 3, characterized in that, Step S5 involves correcting the initial dose contribution matrix, including: A global modulation factor is introduced for each sub-beam j. The global modulation factor is solved by fitting the difference between the Monte Carlo dose contribution matrix and the initial dose contribution matrix using the least squares method. The corrected dose contribution matrix is ​​generated as follows: ； Alternatively, voxels can be divided into several groups according to tissue type, and the correction factor for each sub-beam j within each group can be the average of the correction factors of all voxels within that group. Correction dose contribution matrices can be generated according to tissue type; the tissue types include lung tissue, soft tissue, and bone tissue. Alternatively, a pre-trained machine learning model can be used to predict a correction factor field using the initial dose contribution matrix and local features of the patient image data as input. The initial dose contribution matrix can then be corrected based on the predicted correction factor field to generate a corrected dose contribution matrix.

8. A method for generating radiotherapy plans based on adaptive Monte Carlo calculation according to claim 2 or 3, characterized in that, Step S6, which involves continuing iterative optimization based on the corrected dose contribution matrix until the final convergence condition is met, includes: The current weight vector when the Monte Carlo calculation is triggered As the initial value, the modified dose contribution matrix is ​​used. As the dose calculation engine, iterative optimization continues. In each iteration, the cumulative dose of each voxel v is calculated using the following formula. : ； In the formula, To correct the element values ​​of the dose contribution matrix, Let be the weight of the j-th sub-beam; Calculate the objective optimization function and its gradient based on the cumulative dose, and update the weight vector based on the gradient; When the objective function converges to the preset accuracy Or the number of iterations reaches the preset maximum number of iterations. When the iteration stops, the optimal weight vector is obtained; Optionally, after the final convergence condition is met, step S6 further includes: Independent Monte Carlo simulations were performed based on the obtained optimal weight vectors to obtain the validation dose distribution. The verified dose distribution is compared with the optimized dose distribution calculated based on the modified dose contribution matrix; When the deviation between the two exceeds the preset deviation threshold, the steps of calculating the Monte Carlo dose contribution matrix, matrix correction, and post-correction optimization iteration are repeatedly executed with the current optimal weight vector as the initial value, and the correction is performed twice or more in a loop until the deviation is lower than the preset deviation threshold.

9. A method for generating radiotherapy plans based on adaptive Monte Carlo calculation according to claim 2 or 3, characterized in that, The construction of the initial dose contribution matrix in step S1 includes: Initialize the beam sequence based on the prescribed dosage and target area geometry model; Each beam in the beam sequence is discretized into several sub-beam grids, and each sub-beam is an independent irradiation unit. The total number of sub-beams is N and the total number of voxels is V. The pencil beam analytical dosimetry algorithm is used to calculate the unit dose contribution coefficient of each sub-beam j to each voxel v. ; An initial dose contribution matrix is ​​constructed using the unit dose contribution coefficient of each sub-beam as elements. Each column corresponds to the unit dose contribution distribution of a sub-beam to all voxels, and each row corresponds to the unit dose contribution of a voxel to all sub-beams. The initial dose contribution matrix is ​​stored in the following sparse storage manner: a CSR format sparse index structure is established for the non-zero dose contribution relationship of each sub-beam to the voxel, and only the element values ​​exceeding the preset dose threshold are stored; The cumulative dose of voxel v under any sub-beam weight vector x Calculated using the following formula: 。 10. A method for generating radiotherapy plans based on adaptive Monte Carlo calculation according to claim 1, characterized in that, Step S7, which involves outputting a radiotherapy plan based on the optimal weight vector, includes: Based on the optimal weight vector, the weight distribution of each sub-beam is mapped to the blade motion trajectory of the multi-leaf collimator using a blade sorting algorithm, generating an executable sub-field sequence. Using a Monte Carlo engine independent of step S4, the final dose verification is performed based on the optimal weight vector to obtain the deliverable dose distribution; The deliverable dose distribution is compared with the optimized dose distribution calculated based on the modified dose contribution matrix using dose-volume histogram indices and gamma analysis. When the gamma analysis pass rate meets the preset standard, a treatment plan report is generated; the treatment plan report includes a dose-volume histogram curve, a dose distribution map, a gamma analysis report, and the subfield sequence.

11. A radiotherapy plan generation system based on adaptive Monte Carlo computation, characterized in that, include: Acquisition Module: Used to acquire patient image data, establish geometric models of the target area and organs at risk, set prescription dose and dose constraints, use analytical dose algorithms to calculate the unit dose contribution of each sub-beam to the voxel, and construct the initial dose contribution matrix; The analytical optimization module is used to perform analytical iterative optimization of the dose distribution based on the initial dose contribution matrix, with the weights of each sub-beam as optimization variables, and using an iterative optimization algorithm. In each iteration, the cumulative dose is calculated based on the current weight vector and the initial dose contribution matrix, and the optimized state parameters are calculated based on the cumulative dose. Triggering Decision Module: Used to monitor the convergence trend of the optimization state parameters during the analytical iterative optimization process, and determine whether the analytical optimization has entered a convergence stagnation state based on the convergence trend; when it is determined that it has entered a convergence stagnation state, the Monte Carlo dose calculation is triggered; otherwise, the analytical optimization module is triggered to continue executing the analytical iterative optimization. Monte Carlo calculation module: When Monte Carlo dose calculation is triggered, it uses the Monte Carlo algorithm to calculate the unit dose contribution coefficient of each sub-beam to the voxel and constructs the Monte Carlo dose contribution matrix. Matrix correction module: used to correct the initial dose contribution matrix based on the difference between the initial dose contribution matrix and the Monte Carlo dose contribution matrix, and generate a corrected dose contribution matrix; Continue optimization module: used to continue iterative optimization based on the corrected dose contribution matrix to obtain the optimal weight vector; Output module: Used to output a radiotherapy plan based on the optimal weight vector.

Images