A monte carlo fast dose calculation method

Abstract

The application discloses a kind of monte carlo fast dose calculation methods, load simulation parameter;Get the particle source phase space data of radiotherapy accelerator;Get human anatomy image data, and reconstruct adult three-dimensional matrix of body;Three-dimensional body matrix is gridded;Input monte carlo dose calculation parameter, construct fine monte carlo dose calculation physical model;The application utilizes the structured sparse characteristics of NVIDIAGPU, dynamically compresses input data, so that the effective utilization rate of video memory is improved, and the continuous operation of the computing task is ensured;MIG instance management technology combines multidimensional scheduling and cross-instance collaboration to improve GPU resource utilization, especially in the multi-task parallel scene. The two work together to fully utilize the hardware performance of NVIDIA Ampere and newer generation architecture GPU, improve the efficiency of Monte Carlo dose calculation to a new level, reduce the calculation time, and effectively support the GPU Monte Carlo calculation requirements of particle transport.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of radiotherapy dose calculation, in particular to a Monte Carlo fast dose calculation method. BACKGROUND

[0002] Monte Carlo dose calculation technique is a numerical simulation method relying on random sampling. When dealing with particle transport problems, this technique uses computer-generated pseudo-random numbers to randomly sample relevant parameters in the process of particle physics transport according to a certain probability distribution, thereby simulating the motion trajectory of a single particle in the medium. As a recognized high-precision calculation method in the field of radiotherapy dose calculation, Monte Carlo technique occupies an important position in medical physics and radiation dosimetry due to its accuracy in particle transport simulation and geometric modeling.

[0003] However, although general-purpose Monte Carlo simulation software such as EGS4, EGSnrc, GEANT, MCNP, etc. uses high-precision transport algorithms and reaction cross-section data to simulate the transport process of particles of different energies in various media, the problem of long calculation time is still the main bottleneck in clinical application. Over the years, in order to improve the calculation efficiency and shorten the calculation time, researchers have proposed various fast Monte Carlo calculation models based on different principles, which have significantly improved the calculation speed of the Monte Carlo method. However, to achieve the requirement of near real-time calculation in clinical application, further technical breakthroughs are still needed.

[0004] Since NVIDIA launched the parallel computing architecture CUDA based on GPU in 2006, using GPU for Monte Carlo simulation has gradually become an important direction in dose calculation research. Compared with traditional CPU computing, GPU has significant advantages in cost and energy consumption per computing capacity, which makes it have great application potential in the field of medical physics. During this period, global research teams have developed GPU-based Monte Carlo simulation tools such as gDPM[2], GPUMCD, GMC, SMC, etc.[3-5]. These tools have significantly improved the calculation speed while maintaining the high accuracy of Monte Carlo simulation.

[0005] NVIDIA's Ampere architecture is a revolutionary GPU architecture introduced in 2020, representing a major breakthrough in performance, energy efficiency, and functionality for GPU technology. Based on TSMC's 7nm or Samsung's 8nm process, the architecture significantly improves transistor density and energy efficiency, providing strong hardware support for high-performance computing, artificial intelligence, and graphics rendering. The core innovations of the Ampere architecture include the third-generation Tensor Core and the second-generation RT Core. The third-generation Tensor Core supports multiple data types (such as FP64, TF32, FP16, etc.), significantly accelerating deep learning training and inference tasks, while introducing sparse computing technology to further improve AI computing efficiency by utilizing the sparsity in neural networks. The second-generation RT Core improves real-time ray tracing performance through improved algorithms and hardware acceleration, bringing more realistic visual effects to game development, film rendering, and industrial design.

[0006] The structural sparsity and MIG features of NVIDIA Ampere architecture are important innovations that improve GPU performance and utilization. The structural sparsity feature mainly uses pruning technology to set non-important parameters in neural networks to zero, generating a sparse network. In the NVIDIA Ampere architecture, this sparsity is presented in a 2:4 pattern, meaning that at least two out of every four elements must be zero. This pattern reduces the data footprint and bandwidth of matrix multiplication operations by half, and by skipping the calculation of zero values through sparse Tensor Core, the theoretical computing throughput is doubled. This feature is particularly useful in deep learning inference, as it can significantly accelerate matrix multiplication operations and improve model inference speed while maintaining model accuracy. The MIG (Multi-Instance GPU) feature allows a physical GPU to be securely divided into up to seven independent GPU instances, each with independent memory paths, caches, and computing cores. This division ensures resource isolation for each instance, allowing multiple users or applications to use GPU resources simultaneously without interfering with each other. The MIG feature is particularly suitable for workloads that are not fully saturated, maximizing GPU utilization by running different workloads in parallel. In addition, MIG supports multiple deployment configurations, including bare-metal, virtual machine passthrough, and MIG-based vGPU, providing flexibility for different application scenarios.

[0007] The gDPM, GPUMCD, GMC and other GPU memory-based Monte Carlo programs do not fully utilize the advantages of the Ampere architecture of the new generation of NVIDIA, and the parallelization simulation of the Monte Carlo has further space for mining in the calculation speed. The method designs a set of calculation method adapting to the Ampere and the next generation architecture on the basis of the original Monte Carlo simulation, so as to further improve the calculation speed of the GPU Monte Carlo simulation. SUMMARY

[0008] The purpose of the present application is to provide a Monte Carlo fast dose calculation method which designs a set of calculation method adapting to the NVIDIA Ampere and the next generation architecture graphics card on the basis of the original Monte Carlo simulation, so as to further improve the calculation speed of the GPU Monte Carlo simulation.

[0009] To achieve the above purpose, the present application provides the following technical scheme: comprising the following steps:

[0010] S1: loading simulation parameters;

[0011] S2: obtaining particle source phase space data of a radiotherapy accelerator;

[0012] S3: obtaining human anatomy image data and reconstructing a three-dimensional matrix of the human body;

[0013] S4: meshing the three-dimensional human body matrix;

[0014] S5: inputting Monte Carlo dose calculation parameters and constructing a refined Monte Carlo dose calculation physical model;

[0015] S6: using the structured sparsity of the Ampere architecture to dynamically compress the low-contribution voxels in the input matrix, reduce the memory occupation, data transmission overhead and calculation, and using the "index-value" sparse storage mode for the voxels, assuming that M is the length of the matrix, N is the width of the matrix, and L is the height of the matrix, establishing an index array to record the spatial coordinates (i, j, k), converting the original voxel matrix V into a sparse matrix, and dynamically compressing the low-contribution voxels in the input matrix to reduce the memory occupation, data transmission overhead and calculation amount;

[0016] S7: using the MIG feature of the Ampere architecture, inputting the data of each dose calculation unit using multi-instance MIG calculation, and distributing to different memory, and storing the corresponding particle simulation data in different MIG instance memories in the same way;

[0017] S8: starting the particle transport calculation in different MIG instances, and obtaining the initial energy E and direction V of the particles by sampling;

[0018] S9: After determining the energy E and direction V of the good particles, the collision reaction of the particles in different cross sections and the average free path of the particles are determined by using sampling method combined with the particle cross section data of different materials, the dose deposition on the travel route is calculated, and the result is saved in the corresponding MIG instance memory;

[0019] S10: The normalized grid dose calculation results in each MIG instance are superimposed to obtain the total radiation dose;

[0020] S11: The result is copied to the CPU buffer, and the data in the buffer is saved to the disk.

[0021] Preferably, in S5, the step length is sampled according to the linear attenuation coefficient of the medium, and the step length is randomly sampled by a probability density function, which is closely related to the medium material, particle type and energy, and the value of the corresponding medium is obtained in real time by table lookup method.

[0022] Preferably, in S5, in the interaction type determination link, based on the current energy E of the particle, the atomic number, density and other characteristics of the medium material are combined to call the pre-established cross section database, and the occurrence probability of collision, Compton scattering, photoelectric effect and electron pair production is determined by probability calculation.

[0023] Preferably, in S5, when calculating the energy deposition, the energy deposition amount ΔE is calculated by using the corresponding physical formula according to the specific interaction type, in the photoelectric effect, the energy deposition is related to the electron binding energy of the atom and the photon energy, and the residual energy of the particle is updated by formula (Bi is the electron binding energy) calculation. .

[0024] Preferably, in S6, the voxel contribution degree calculation model is constructed, for the voxels in the three-dimensional dose calculation space, the dose contribution degree C of the voxel is determined by the energy deposition accumulation of the particle transport, wherein N is the number of particles passing through the voxel, ΔEn is the energy deposition of particle n, δ is a position determination function, a contribution degree threshold Dth is set, voxels below the threshold are determined as low contribution voxels, wherein is a coefficient between 0 and 1, voxels below the contribution degree threshold are removed, voxels above the contribution degree threshold are used, voxels with less influence on the particle collision process are considered, and the memory occupation and calculation amount are reduced.

[0025] Preferably, the S7 detects the graphic card hardware resource, including the SM unit quantity, the graphic memory capacity, the graphic memory bandwidth, the register file size, checks whether the graphic card supports the MIG feature, loads the customized computing kernel program in the instance initialization stage, initializes the graphic memory space, establishes the thread scheduling queue inside the instance, ensures the efficient execution of the computing task, configures different graphic memory occupation units according to the size of the input data graphic memory occupation, and reserves part of the shared resources for data interaction.

[0026] Preferably, after the particle simulation in S10 is completed, a computing result is obtained in each MIG instance, the results of the dose calculation in the computing area are superimposed, and the total radiation dose distribution is obtained.

[0027] In summary, due to the adoption of the above-mentioned technology, the beneficial effects of the present application are:

[0028] Compared with the prior art, the embodiment of the present application provides a GPU parallel Monte Carlo fast dose calculation method based on the NVIDIA architecture, which optimizes the transport process of the Monte Carlo simulation particles through two main methods. The present application has the following advantages: the structured sparse characteristics of the NVIDIA graphic card are utilized to dynamically compress the input data, so that the effective utilization rate of the graphic memory is improved, and the continuous operation of the computing task is ensured; the MIG instance management technology is combined with multi-dimensional scheduling and cross-instance cooperation, so that the GPU resource utilization rate is improved, especially in the multi-task parallel scene. The two work together to fully exert the hardware performance of the NVIDIA Ampere and the newer generation of architecture graphic card, improve the Monte Carlo dose calculation efficiency to a new level, reduce the calculation time, and effectively support the GPU Monte Carlo computing demand of particle transport. BRIEF DESCRIPTION OF DRAWINGS

[0029] The accompanying drawings, which form a part of this application, are included to provide a further understanding of the application and are incorporated in and constitute a part of this application. The illustrative embodiments of the application are described and explained with additional specificity and detail through the use of the accompanying drawings in which:

[0030] Figure 1 The accompanying drawings, which form a part of this application, are included to provide a further understanding of the application and are incorporated in and constitute a part of this application. The illustrative embodiments of the application are described and explained with additional specificity and detail through the use of the accompanying drawings in which: DETAILED DESCRIPTION

[0031] 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 a part of the embodiments of the present invention, not all of them. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention. Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to represent selected embodiments of the invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0032] In the description of this invention, it should be understood that the terms indicating orientation or positional relationship are based on the orientation or positional relationship shown in the drawings and are only for the convenience of describing the invention and simplifying the description, and are not intended to 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 the invention.

[0033] In this invention, unless otherwise explicitly specified and limited, the terms "installation," "connection," "linking," and "fixing," etc., should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral part; they can refer to a mechanical connection or an electrical connection; they can refer to a direct connection or an indirect connection through an intermediate medium; they can refer to the internal communication of two components or the interaction between two components. Those skilled in the art can understand the specific meaning of the above terms in this invention based on the specific context of the specification.

[0034] This invention provides a Monte Carlo rapid dose calculation method, a parallel Monte Carlo rapid dose calculation method based on NVIDIA Ampere graphics card, comprising the following steps.

[0035] S1: Load simulation parameters

[0036] S2: Acquire particle source phase space data for radiotherapy accelerators

[0037] S3: Acquire human anatomical structure image data and reconstruct a three-dimensional matrix of the human body.

[0038] S4: Mesh the 3D human body matrix

[0039] S5: Input Monte Carlo dose calculation parameters and construct a refined physical model for Monte Carlo dose calculation:

[0040] For particle motion step size sampling, based on the linear attenuation coefficient of the medium Step length by probability density function Random sampling of the step length, wherein The value is closely related to the medium material, particle type and energy, and is obtained in real time by table lookup method

[0041] In the interaction type determination link, based on the current energy E of the particle, the atomic number, density and other characteristics of the medium material are combined to call the pre-established cross-section database, and the occurrence probability of collision, Compton scattering, photoelectric effect, electron pair production and other interactions is determined by probability calculation;

[0042] In the energy deposition calculation, the corresponding physical formula is used to calculate the energy deposition amount ΔE according to the specific interaction type. For example, in the photoelectric effect, the energy deposition is related to the electron binding energy of the atom and the photon energy, which is calculated by the formula (Bi is the electron binding energy) to update the remaining energy of the particle .

[0043] S6: Using the structured sparsity of Ampere architecture, low-contribution voxels in the input matrix are dynamically compressed to reduce memory usage, data transmission overhead and calculation. The "index-value" sparse storage mode is used for voxels, assuming that M is the length of the matrix, N is the width of the matrix, and L is the height of the matrix. The index array records its spatial coordinates (i, j, k), and the original voxel matrix V ∈ R (M × N × L) is converted into a sparse matrix. Dynamic compression is performed on low-contribution voxels in the input matrix, effectively reducing memory usage, data transmission overhead and calculation.

[0044] First, a voxel contribution calculation model is constructed. For a voxel in the three-dimensional dose calculation space, its dose contribution C is determined by the energy deposition accumulation of particle transport:

[0045]

[0046] Where N is the number of particles passing through the voxel, ΔEn is the energy deposition of particle n, and δ is the position determination function (δ=1 when the particle is located at (i,j,k), otherwise 0). Set the contribution threshold Dth, and the voxels below the threshold are determined as low-contribution voxels.

[0047]

[0048] Where is a coefficient between 0 and 1. Set the contribution threshold below to be removed, and only use voxels above the contribution threshold, so as to only consider voxels that have less impact on the particle collision process, reducing memory usage and calculation. ​

[0049] S7: Using the MIG feature of the Ampere architecture, input the data of each dose calculation unit into the calculation of multiple instances of MIG, and distribute them to different video memories, and store the corresponding particle simulation data in the same way in the different MIG instance video memories

[0050] First, detect the hardware resources of the graphics card, including the number of SM units, the capacity of the video memory, the bandwidth of the video memory, the size of the register file, etc. Check whether the graphics card supports the MIG feature. In the instance initialization stage, load the customized calculation kernel program, initialize the video memory space, establish the thread scheduling queue inside the instance, and ensure efficient execution of the calculation task.

[0051] According to the size of the input data video memory occupation, configure different video memory occupation units, and reserve part of the shared resources for data interaction. Therefore, the MIG instance division needs to meet the video memory capacity constraint condition, and reserve space for dynamic adjustment.

[0052] We first set up several core variables:

[0053] Input data video memory Vinput: input data, medium parameters, etc. Static data occupation

[0054]

[0055] Particle trajectory dynamic video memory Vdynamic: real-time generated particle state (position / energy / history)

[0056]

[0057] Compute reserved space Vcompute: intermediate result cache, thread stack, burst memory growth buffer, store intermediate data such as scattering cross section temporary table and energy deposition accumulation array

[0058]

[0059] Dynamic adjustment factor: safety factor α (1.2~1.5), to cope with local video memory peak caused by particle clustering, based on historical task variance dynamic adjustment

[0060] Memory occupation constraint of a single MIG

[0061]

[0062] The number of MIGs is:

[0063]

[0064] Downward rounding, to ensure that the number of MIG instances is an integer and can meet the needs of all capacities, The maximum number of instances supported by the graphics card.

[0065] In a single MIG instance, the number of SM units is allocated according to the set number of particles.

[0066] S8: In different MIG instances, the transport calculation of particles is started, and the initial energy E and direction V of the particles are obtained by sampling.

[0067] S9: After determining the energy E and direction V of the particles, the sampling method is used to combine the particle cross-section data of different materials to determine the collision reaction of particles in different cross-sections and the average free path of particle travel, calculate the dose deposition on the travel route, and save the results to the corresponding MIG instance memory.

[0068] S10: The normalized grid dose calculation results in each MIG instance are superimposed to obtain the total radiation dose.

[0069] After the particle simulation is completed, a calculated result is obtained in each MIG instance, and finally the results of multiple dose calculations in the calculation region are superimposed to obtain the total radiation dose distribution

[0070]

[0071] where n is the number of MIGs, Wn is the weight of the nth MIG result, is the result obtained by the nth MIG instance, is the total calculation result.

[0072] S11: The results are copied to the CPU buffer, and then the data in the buffer are saved to the disk.

[0073] Embodiment:

[0074] On the basis of existing GPU Monte Carlo calculation, the transport process of Monte Carlo simulation particles is optimized by two main methods to improve the calculation efficiency of the GPU Monte Carlo method and reduce the time required for calculation. One is dynamic compression optimization for low-contribution voxels in input data, which only considers voxels that have a greater impact on particles in the calculation of particle transport calculation, thereby reducing the memory occupation and calculation of invalid voxels. The second method is to use the MIG technology of the NVIDIA Ampere architecture to perform fine-grained configuration of instances, divide the original Monte Carlo program running on a single GPU into multiple independent entities, and perform calculation respectively, thereby improving the calculation speed.

[0075] The method specifically includes the following steps (as shown in Figure 1 )

[0076] Load the required data, including simulation data, particle cross-section data, etc. The loaded simulation data can be human CT data, measured phantom data, or other solid data of different materials. The phase space data loaded can also be the spatial distribution formula of particles, including energy spectrum, flux, etc., or specific structural materials and shapes.

[0077] Input the Monte Carlo dose calculation parameters, build a refined Monte Carlo dose calculation physical model, and build a highly refined particle transport model based on the basic principles of Monte Carlo method. In the model building, the key physical processes of particle motion in the medium are accurately defined:

[0078] In Monte Carlo simulation, particle motion step length is realized by exponential distribution sampling: according to the linear attenuation coefficient of the medium (get the μ value of the current medium by table lookup), generate a random step length that satisfies the probability density function, simulate the random collision behavior of particles in the medium. This process uses hardware-accelerated inverse transform sampling instructions to ensure that the calculation time per million particle steps is <5ms, providing the basis for the spatial and temporal distribution of dose deposition.

[0079] Energy deposition calculation is based on differentiated processing of reaction types: in photoelectric effect, energy deposition (Bi is the electron binding energy), and Compton scattering calculates energy loss according to the scattering angle formula to update the remaining energy of the particle . This process uses Ampere's Tensor Core to accelerate floating-point operations, ensuring that the single-particle energy update delay is <20ns, supporting real-time tracking of millions of particles.

[0080] In this embodiment, the photon reaction is taken as an example, in addition to which, the reactions of protons, carbon ions and other related particles are similar in principle. Those skilled in the art should know the processes of collision and reaction of different particle types in different materials, and this embodiment does not limit specific reaction particles and transport scenarios.

[0081] Based on the structure sparsification characteristics of NVIDIA Ampere architecture, Monte Carlo simulation is commonly used in dose calculation and particle transport simulation in radiotherapy. These processes involve a large number of matrix operations. Using the structure sparsification compression technology of NVIDIA graphics card, these matrices can be pruned and compressed to reduce the amount of calculation and storage requirements, thereby improving the efficiency and speed of simulation.

[0082] The embodiment is directed to dynamic compression optimization of low contribution voxels in dose calculation. Ampere architecture supports structured sparse computation, based on which, an "index-value" sparse storage mode is used for voxels, assuming that M is the length of the matrix, N is the width of the matrix, and L is the height of the matrix, an index array is established to record the spatial coordinates (i, j, k), and the original voxel matrix V is converted into a sparse matrix, and dynamic compression is implemented for low contribution voxels in the input matrix, thereby effectively reducing memory occupation, data transmission overhead and calculation amount. The voxels use an "index-value" sparse storage mode, and the sparse storage mode can be a one-dimensional index value such as D(i)=Value, or a three-dimensional index value D(i,j,k)=Value, and different memory spaces are allocated according to different settings.

[0083] In an example embodiment of the application, the determination of low contribution voxels is determined by the formula in the summary, and the determination of low contribution voxels can also be obtained by calculating one or a combination of physical factors and biomedical factors based on the method of the user's region of interest; wherein the physical factors are physical conditions reflecting the material composition of the patient or phantom and the irradiation; the material composition of the patient or phantom includes the density, CT value, mass number and atomic number of the phantom; the irradiation physical conditions include: field distribution, source distribution. The biomedical factors include: organ tissue irradiation threshold, biological sensitivity damage probability, etc.

[0084] In addition, different application scenarios have obvious differences in corresponding input phantom data, and the voxel values of different phantoms have obvious differences in particle collision response, and different threshold settings can be selected to remove voxels. Taking CT as an example, if the CT gray value is less than 1000, we can understand it as air, which has very small collision, blocking and absorption of photons or protons, and the voxels with CT gray value less than 1000HU can be removed, and only the voxels with CT gray value greater than or equal to 1000HU are retained.

[0085] Using the MIG feature of the Ampere architecture, the data input of each dose calculation unit is calculated using multiple instance MIG, and is distributed to different video memories, and the corresponding particle simulation data is stored in the same way in different MIG instance video memories.

[0086] First, detect the hardware resources of the graphics card, including the number of SM units, the capacity of the video memory, the bandwidth of the video memory, the size of the register file, whether the MIG feature is supported, etc. In the instance initialization stage, load the customized computing kernel program, initialize the video memory space, establish the thread scheduling queue inside the instance, and ensure efficient execution of the computing task.

[0087] According to the size of the input data memory occupation, different memory occupation units are configured, and part of the shared resources are reserved for data interaction. Therefore, the MIG instance division needs to meet the memory capacity constraint condition, and at the same time, a dynamic adjustment space is reserved. The related calculation formula is referred to in the part of the invention.

[0088] In different MIG instances, the transport calculation of the starting particles is started. First, the initial energy E and direction V of the particles are obtained by sampling. After the energy E and direction V of the particles are determined, the sampling method is used to combine the particle cross-section data of different materials to judge the collision reaction of the particles in different cross-sections and the average free path of the particles. The dose deposition on the travel route is calculated, and the result is saved to the memory of the corresponding MIG instance. The normalized grid dose calculation results in each MIG instance are superimposed to obtain the total radiation dose. After the particle simulation is completed, a calculation result is obtained in each MIG instance. Finally, the results of multiple dose calculations in the calculation region are superimposed to obtain the total radiation dose distribution. Finally, the results are copied to the CPU buffer, and then the data in the buffer is saved to the disk.

[0089] In the Monte Carlo simulation, the MIG (Multi-Instance GPU) technology of the NVIDIA Ampere architecture is used to split the calculation task, which can realize hardware-level resource isolation and parallel optimization: by dividing a single card into independent instances (such as 10GB memory / 10SM or 40GB memory / 40SM), each instance exclusively occupies memory and calculation units, avoiding resource competition between multiple tasks and improving stability; at the same time, multiple cases can be processed in parallel (such as 3 MIG instances calculating 3 fields simultaneously), and cross-instance particle data interaction is realized through NVLink high-speed communication, with a throughput increase of 3-5 times; dynamic memory slicing (such as each instance storing only 20% of the dose matrix) combined with structure sparsification compression reduces memory occupation by 60%, supporting real-time calculation of complex cases; the efficiency is 2 to 3 times higher than that of the traditional single-instance mode.

[0090] In addition, in 2022, NVIDIA launched a new generation of acceleration computing platform Hopper architecture, which is designed for large-scale AI, HPC and supercomputer workloads, also has structure sparsification and MIG characteristics. Therefore, the present application is also applicable to Hopper architecture graphics cards and all graphics cards supporting these two characteristics released by NVIDIA in the future.

[0091] Compared with the prior art, the application has the following advantages: the input data is dynamically compressed by using the structural sparsification characteristics of the NVIDIA graphics card, the effective utilization rate of the display memory is improved, and the continuous operation of the calculation task is ensured; the MIG instance management technology is combined with multi-dimensional scheduling and cross-instance cooperation, so that the utilization rate of GPU resources is improved, and the advantages are obvious especially in the multi-task parallel scene. The two work together to fully exert the hardware performance of the NVIDIA Ampere and newer generation architecture graphics card, and improve the efficiency of the Monte Carlo dose calculation to a new level, and effectively support the Monte Carlo calculation demand of particle transport.

Claims

1. A Monte Carlo fast dose calculation method, characterized by: The method comprises the following steps: S1: loading simulation parameters; S2: obtaining particle source phase space data of a radiotherapy accelerator; S3: obtaining human anatomical structure image data and reconstructing a human three-dimensional matrix; S4: meshing the human three-dimensional matrix; S5: inputting Monte Carlo dose calculation parameters and constructing a refined Monte Carlo dose calculation physical model; S6: using the structured sparsity of the Ampere architecture to dynamically compress low-contribution voxels in the input matrix, reduce memory occupation, data transmission overhead and calculation, and using the "index-value" sparse storage mode for voxels, assuming that M is the length of the matrix, N is the width of the matrix, and L is the height of the matrix, an index array records the spatial coordinates (i, j, k), and the original voxel matrix V is converted into a sparse matrix, and low-contribution voxels in the input matrix are dynamically compressed; S7: using the MIG feature of the Ampere architecture, using a multi-instance MIG calculation method, inputting data of each dose calculation unit for processing, distributing the data to different memories, and storing corresponding particle simulation data in different MIG instance memories in the same way; S8: in different MIG instances, starting particle transport calculation, and obtaining initial energy E and direction V of the particles through sampling; S9: after the energy E and direction V of the particles are determined, the sampling method is used, the particle cross-section data of different materials are combined, the collision reaction of the particles in different cross-sections and the average free path of the particles are judged, the dose deposition on the travel route is calculated, and the result is saved in the memory of the corresponding MIG instance; S10: superimposing the normalized grid dose calculation results in each MIG instance to obtain the total radiation dose; S11: copying the total radiation dose to the CPU buffer, and saving the data in the buffer to the disk.

2. The Monte Carlo fast dose calculation method of claim 1, wherein: In S5, the step length is sampled according to the linear attenuation coefficient of the medium, and the step length is randomly sampled through the probability density function, the sampling value is closely related to the medium material, the particle type and the energy, and the value of the corresponding medium is obtained in real time through the table lookup method.

3. The Monte Carlo fast dose calculation method of claim 2, wherein: In S5, in the interaction type judgment link, based on the current energy E of the particle, the atomic number and density characteristics of the medium material are combined, the pre-established cross-section database is called, and the occurrence probability of collision, Compton scattering, photoelectric effect and electron pair production interaction is determined through probability calculation.

4. The Monte Carlo fast dose calculation method of claim 3, wherein: In S5, when the energy deposition is calculated, the energy deposition amount ΔE is calculated according to the specific interaction type, in the photoelectric effect, the energy deposition is related to the electron binding energy of the atom and the photon energy, ΔE=E-Bi is calculated, Bi is the electron binding energy, and the remaining energy of the particle is updated as E'=E-ΔE.

5. The Monte Carlo fast dose calculation method of claim 1, wherein: The S6 voxel contribution degree calculation model is used to calculate the dose contribution degree C of a voxel in a three-dimensional dose calculation space by energy deposition accumulation of particle transport, wherein N is the number of particles passing through the voxel, ΔEn is the energy deposition of particle n, and δ is a position determination function. A contribution degree threshold Dth is set, and voxels below the threshold are determined as low contribution voxels, wherein is a coefficient between 0 and 1. The voxels below the contribution degree threshold are removed, the voxels above the contribution degree threshold are used, the voxels affected by particle collision are considered, and the memory occupation and calculation amount are reduced.

6. The Monte Carlo fast dose calculation method of claim 1, wherein: In S7, the hardware resources of the display card are detected, including the number of SM units, the memory capacity, the memory bandwidth, the register file size, and whether the display card supports the MIG feature. In the instance initialization stage, a customized calculation kernel program is loaded, the memory space is initialized, the thread scheduling queue in the instance is established, the efficient execution of the calculation task is ensured, different memory occupation units are configured according to the size of the input data memory occupation, and part of the shared resources are reserved for data interaction. The MIG instance division needs to meet the memory capacity constraint condition, and a dynamic adjustment space is reserved.

7. The Monte Carlo fast dose calculation method of claim 1, wherein: In S10, after the particle simulation is completed, a calculation result is obtained in each MIG instance. The results of multiple dose calculations in the calculation region are superimposed to obtain the total radiation dose distribution.

Images