BNCT dose calculation method, electronic equipment and computer readable storage medium

Abstract

The invention discloses a BNCT dose calculation method, which comprises the following steps: S1, constructing a voxel model, initializing one or more particles in a simulation system, and respectively acquiring the initial position, energy and movement direction of each particle; s2, adaptively selecting a Monte Carlo particle tracking algorithm and a flux counting mode based on the region where the particles are located and the energy of the particles; s3, repeating the step S2 until the particles leave the simulation system or are killed when meeting a truncation condition; and S4, cumulatively calculating dose distribution. Through the BNCT dose calculation method, the Monte Carlo particle transport simulation time is greatly shortened, the BNCT dose calculation efficiency is optimized, and the method is suitable for rapid and high-precision radiation dose simulation under a BNCT complex tissue structure.

Description

Technical Field

[0001] This invention relates to the fields of nuclear technology applications and medical physics, specifically to a BNCT dose calculation method, electronic device, and computer-readable storage medium. Background Technology

[0002] Boron neutron capture therapy (BNCT) is an advanced, cell-level targeted radiotherapy technique. The formulation of BNCT treatment plans relies on precise dose calculations. Because BNCT dose calculations involve simulating neutron / photon coupled particle transport within the complex and inhomogeneous structures of the human body, Monte Carlo particle transport procedures are typically used as the core of the dose calculation. In Monte Carlo particle transport simulations, one of the core algorithms is the particle tracking algorithm, which determines the random walk process and state of particles within the simulated system. Currently, the mainstream particle tracking algorithms include Ray Tracking and Delta Tracking.

[0003] In existing technologies, Ray Tracking is commonly used for particle tracking, and dose statistics are calculated using track length counting. Each particle contributes to the count for every voxel grid it passes through. However, within the voxel model, frequent boundary crossing calculations result in low efficiency and long dose calculation times for this particle tracking algorithm. Delta Tracking treats the entire simulation system as a single medium for particle tracking, which greatly reduces the calculation time for particle boundary crossings. However, it can only use collision counting, leading to larger statistical errors in dose. It requires simulating more particles to ensure the accuracy of the results, and the calculation time is not reduced, making it difficult to meet the needs of clinical BNCT. Summary of the Invention

[0004] The purpose of this invention is to provide a BNCT dose calculation method to solve the problem that traditional particle tracking algorithms are not very efficient when dealing with complex voxel model geometry in existing BNCT dose calculations.

[0005] To solve the above-mentioned technical problems, the present invention adopts the following technical solution:

[0006] In a first aspect, the present invention provides a BNCT dose calculation method, comprising the following steps:

[0007] S1. Construct a voxel model and initialize one or more particles, obtaining the initial position, energy, and direction of motion of each particle;

[0008] S2. Based on the region where the particle is located and the energy of the particle, adaptively select the Monte Carlo particle tracking algorithm and flux counting method to count the dose;

[0009] S3. Repeat step S2 until the particle leaves the simulation system or is killed by meeting the cutoff condition;

[0010] S4. Cumulative calculation of dose distribution.

[0011] Preferably, in step S1, the voxel model is a digital model, obtained by reading information from medical imaging files such as CT, MRI, and PET / CT of the case and converting it into a format that can be recognized by the Monte Carlo particle transport program.

[0012] Preferably, in step S2, the implementation of the adaptive selection of the Monte Carlo particle tracking algorithm and the flux counting method includes:

[0013] S21. Determine the region of each particle;

[0014] S22. If the particle is located in the air region, Delta Tracking is used for particle tracking; if the particle is located inside the human body region, proceed to step S23.

[0015] S23. If the energy of the particle is higher than the preset energy threshold E th If the particle tracking is positive, Ray Tracking is used; otherwise, Delta Tracking is used.

[0016] Preferably, in step S3, the cutoff condition is any one of particle energy cutoff, particle weight cutoff, or time cutoff.

[0017] Preferably, in step S4, the step of accumulating the dose distribution includes: calculating the photon equivalent dose of the voxel grid in the voxel model in each voxel grid, and accumulating it into the dose distribution of the voxel model.

[0018] Preferably, the Delta Tracking includes: a cross-sectional value Σ based on a cross-sectional value greater than or equal to the maximum total macroscopic cross-section of all materials in the simulation system. maj To sample the transport length of the particle, and after the particle collides, based on the macroscopic total cross-section of the material at the particle's current location and Σ maj The ratio is used to determine whether the collision is a real physical collision or a virtual collision in a probabilistic way.

[0019] Preferably, the Ray Tracking includes: directly sampling the distance from the particle to the next physical collision point based on the macroscopic total cross section of the material corresponding to the voxel grid where the particle is currently located, calculating the distance from the particle to the next geometric boundary along the current direction, and taking the smaller of the two values ​​as the next track length of the particle.

[0020] In a second aspect, the present invention provides an electronic device including a processor coupled to a memory; the memory is used to store a computer program, and the processor is used to execute the computer program stored in the memory to implement the BNCT dose calculation method as described above.

[0021] Thirdly, the present invention also provides a computer-readable storage medium for storing a computer program to implement the BNCT dose calculation method as described above.

[0022] The present invention has the following gain effects:

[0023] This invention provides a BNCT dose calculation method, electronic device, and computer-readable medium. By determining the region where the particle is located and the relative magnitude of the particle energy and a preset energy threshold, it adaptively switches between two particle tracking algorithms: Ray Tracking and Delta Tracking. When Ray Tracking is selected, track length counting is used, and when Delta Tracking is selected, collision counting is used. This reduces the simulation time of particles in uninterested areas. When the energy is higher than the preset energy threshold, Ray Tracking ensures the efficiency of dose counting. When the energy is lower than the preset energy threshold, it greatly reduces the time for particle penetration calculation, thereby increasing the overall efficiency of BNCT dose calculation. Attached Figure Description

[0024] To make the objectives, technical solutions, and advantages of the invention clearer, the invention will now be described in further detail with reference to the accompanying drawings, wherein:

[0025] Figure 1 This is a flowchart of the coupled Monte Carlo particle tracking algorithm of the present invention.

[0026] Figure 2 This is a comparison chart of Ray Tracking values ​​along the beam center axis in the voxel model and the FOM values ​​of the method proposed in this invention. Detailed Implementation

[0027] The technical solutions of some 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 disclosed in the present invention, and not all of them. All other embodiments obtained by those skilled in the art based on the embodiments disclosed in the present invention are within the scope of protection of the present invention. It should be noted that in the drawings with reference numerals, the same reference numerals and letters represent similar parts. Once a part is defined in a drawing, it will not be defined and explained again in subsequent drawings.

[0028] Before providing a further detailed description of the present invention, the nouns and terms used in the embodiments are explained, and the nouns and terms used in the embodiments of the present invention are subject to the following interpretations:

[0029] (1) Particle tracking algorithm: In the Monte Carlo particle transport procedure, particles move randomly in phase space, and their state S is characterized by parameters such as position r, angle Ω, energy E, and weight W. The particle is determined by its state S. m =(r m ,Ω m E m W m (m=0,1,2,…,M-1) After moving a distance L, a new state S is reached. m+1 =(r m+1 ,Ω m+1 E m+1 W m+1 The process of (m = 1, 2, 3, ..., M) is called the particle tracking process, and the algorithm that implements this process is called the particle tracking algorithm.

[0030] (2) Surface calculation: Calculate the shortest distance from the particle along the current direction of motion to the boundary of the geometry.

[0031] (3) Ray Tracking: The particle transport length is sampled using the macroscopic total cross section at the particle's location, and a cross-surface calculation is performed. The minimum value of the transport length and the closest distance from the particle along the current direction to the boundary of the geometry is selected as the length of this motion, and the position of the particle in the next state is determined.

[0032] (4) Delta Tracking: The entire system is treated as a single material, and the total macroscopic cross-section of this material is made greater than or equal to the total macroscopic cross-section of other materials in the system. The particle transport length is sampled using this total macroscopic cross-section to determine the next collision point. After one transport step, a rejection sample is performed to determine whether the collision is a real or fake collision.

[0033] (5) Mean free path: The reciprocal of the total cross section of the geometric body containing the particle is the mean free path.

[0034] (6) Track length counting: The trajectory of a particle in a geometric body is called the track, and the track length of a particle in the geometric body is the particle flux of this transport step.

[0035] (7) Collision count: When a particle collides with a geometry, it contributes to the particle flux within that geometry.

[0036] (8) Medical imaging: The medical imaging refers to one or more of CT, MRI, PET / CT.

[0037] (9) Voxel model: Based on the geometric, image grayscale, delineation and other information contained in the medical image, a model is constructed consisting of several rectangular voxel meshes of the same size, including the voxel mesh side length, the material of each voxel mesh, the nuclide composition of each material, the material density and other model parameters.

[0038] (10) Photon equivalent dose: Different dose components produce different biological effects on the human body. In order to evaluate its biological effect, the calculated absorbed dose is multiplied by the corresponding relative biological effect / complex biological effect factor and summed to obtain the photon equivalent dose.

[0039] (11) FOM: Figure of Merit, which is equal to the reciprocal of the product of the square of the relative statistical error and the computation time. Its value represents the efficiency of the program computation.

[0040] This invention can be applied to BNCT Monte Carlo dose calculation, solving the problem that traditional particle tracking algorithms are not very efficient when dealing with complex voxel model geometry in existing BNCT dose calculations.

[0041] The technical concept of this invention is as follows: based on the region where the particle is located in the voxel model and the particle energy, an appropriate particle tracking algorithm is automatically selected. When the particle is in the air region, Delta Tracking is used to accelerate particle tracking, thereby reducing the calculation of cross-area in the air. When the particle enters the human body, if the particle energy is higher than the set energy threshold and the mean free path of the particle is longer, Ray Tracking combined with track length counting is used to ensure the efficiency of dose counting. When the energy is lower than the set threshold, Delta Tracking is used to reduce a large amount of cross-area calculation, thereby reducing the particle simulation time, improving particle tracking efficiency, and thus improving the overall efficiency of BNCT dose calculation.

[0042] The BNCT dose calculation method, electronic device, and computer-readable storage medium disclosed in this invention have the following technical effects: achieving rapid and accurate dose calculation when formulating a BNCT treatment plan.

[0043] To further illustrate the method of the present invention, the following application examples are disclosed.

[0044] In some embodiments, a BNCT dose calculation method is provided, including the following steps:

[0045] S1. Construct a voxel model and initialize one or more particles in the simulation system, and obtain the initial position, energy and direction of motion of each particle respectively;

[0046] S2. Based on the region where the particle is located and the energy of the particle, adaptively select the Monte Carlo particle tracking algorithm and flux counting method to count the dose;

[0047] S3. Repeat step S2 until the particle leaves the simulation system or is killed by meeting the cutoff condition;

[0048] S4. Cumulative calculation of dose distribution.

[0049] In some embodiments, in step S1, the voxel model is a digital model, obtained by reading information from medical image files of a case and converting it into a format recognizable by a Monte Carlo particle transport program. This approach can retain a large amount of grayscale and delineation information in medical images, achieving high-precision voxel model modeling.

[0050] In some embodiments, the implementation of the adaptive selection of the Monte Carlo particle tracking algorithm and the flux counting method in step S2 includes:

[0051] S21. Determine the region of each particle;

[0052] S22. If the particle is located in the air region, Delta Tracking is used for particle tracking; if the particle is located inside the human body region, proceed to step S23.

[0053] S23. If the energy of the particle is higher than the preset energy threshold E th If the particle tracking is positive, Ray Tracking is used; otherwise, Delta Tracking is used.

[0054] This scheme allows for the selection of a suitable particle tracking algorithm based on the physical properties of the particles and the physical properties of the spatial region, enabling efficient BNCT dose calculation.

[0055] In some embodiments, in step S3, the cutoff condition is any one of particle energy cutoff, particle weight cutoff, or time cutoff. This scheme allows the cutoff condition to terminate particle simulation in a timely manner, preventing the program from consuming computational resources on particles with low or no contribution. Specifically, particle energy cutoff refers to the situation where the particle energy is lower than a preset cutoff energy E. min When the particle's weight is lower than the preset cutoff weight, a roulette wheel is played. If the roulette wheel fails, the simulation of the current particle is terminated. When the particle's running time is higher than the preset cutoff time, the simulation of the current particle is terminated.

[0056] In some embodiments, step S4, the step of accumulating the dose distribution, includes: calculating the photon equivalent dose of each voxel grid in the voxel model and accumulating it into the dose distribution of the voxel model. Through this scheme, the photon equivalent dose of each voxel grid and the photon equivalent dose distribution of the entire voxel model can be obtained.

[0057] In some embodiments, Delta Tracking includes: a cross-sectional value Σ based on a cross-sectional value greater than or equal to the maximum total macroscopic cross-section of all materials within the simulation system. maj The transport length of the particle is sampled, and after the particle collides, the macroscopic total cross section of the material corresponding to the voxel mesh where the particle is currently located is calculated based on Σ. maj The ratio is used to determine whether the collision is a real physical collision or a virtual collision in a probabilistic way. Through this scheme, the simulation system is regarded as a single material medium for particle transport simulation, avoiding a large amount of cross-surface calculation and saving a lot of transport simulation time.

[0058] In some embodiments, the Ray Tracking includes: directly sampling the distance from the particle to the next physical collision point based on the macroscopic total cross section of the material corresponding to the current voxel grid of the particle, and calculating the distance from the particle to the next geometric boundary along the current direction, taking the smaller of the two values ​​as the next track length of the particle. With this scheme, even if the particle does not collide in the voxel grid, as long as it passes through the voxel grid, flux counting can be generated. When the particle crosses multiple voxel grids before colliding, the counting efficiency is high.

[0059] In some embodiments, please refer to Figure 1 This paper provides a specific application example of a BNCT dose calculation method. The working principle of the technical solution in this embodiment is as follows: based on the region where the particle is located and the energy of the particle, an appropriate particle tracking algorithm and dose counting method are adaptively switched, thereby improving the overall dose calculation efficiency of BNCT. When the particle is located in the air region, Delta Tracking is switched to reduce the calculation of cross-surface in the air and reduce simulation time. When the particle enters the human body, it is necessary to determine the relative magnitude of the particle energy and the preset energy threshold. If the energy of the particle is higher than the preset energy threshold, Ray Tracking is switched to use track length counting, thereby ensuring the efficiency of dose counting when the mean free path of the particle is long. When the energy of the particle is lower than the preset energy threshold, Delta Tracking is switched to use collision counting, which greatly reduces the time of frequent cross-surface calculation of particles in complex voxel models.

[0060] The working process of this solution is as follows:

[0061] 1) Construct a voxel model based on the medical images. The voxel model has a voxel size of x×y×z and contains P types of materials and Q types of nuclides.

[0062] 2) Sample the source particle, including the particle's exit position r, angle Ω, and energy E, to determine the initial state S0(r0,Ω0,E0,W0) of the source particle.

[0063] 3) Determine the type of region the particle is in based on its position r;

[0064] 4) If the particle is located in the air region, switch to the Delta Tracking particle tracking algorithm, calculate the maximum macroscopic total cross-section of material P based on the particle energy E, and assign this value to Σ. maj The entire system is considered as a single material at its macroscopic cross-section. Using Σ maj Sampling of transport length is performed to determine the position r' of the particle's next collision, based on the macroscopic total cross-section of the material of the geometry containing r' and Σ. maj Perform a rejection sample. If it is a true collision, confirm the particle state after the collision and use collision counting to count the dose, then proceed to step 7.

[0065] 5) If the particle is located inside the human body, further determine the particle's current energy E and the preset energy threshold E. th The relative size, if E>E th If the particle energy E is used, the Ray Tracking particle tracking algorithm is switched to. The macroscopic total cross-sectional sampling transport length L of the geometry at the current position r is calculated based on the particle energy E. The minimum distance d from r to the current geometry boundary is calculated. The smaller value of L and d is taken as the track length of the transported particle in the geometry. The particle position is updated and the dose is counted using track length counting. Then, proceed to step 7.

[0066] 6) If E <E th Then switch to the Delta Tracking particle tracking algorithm and use collision counting to count the dose. The algorithm execution process and counting method are the same as in step 4.

[0067] 7) If the particle's history ends, i.e. the particle leaves the simulation system, or the particle meets the truncation conditions such as energy truncation, weight truncation, or time truncation, then stop the simulation of that particle; otherwise, return to step 2 to re-determine the region where the particle is located and perform subsequent simulations.

[0068] 8) If the simulation of all source particles is completed, terminate the program; otherwise, return to step 1 and restart the sampling of source particles.

[0069] The method was tested on a voxel model with a resolution of 512×512×140, where each voxel grid measures 1.2832mm×1.2832mm×3.0000mm. The voxel model contains 43 materials and 16 nuclides. Dose calculations were performed using both Ray Tracking and the above method. The doses and their statistical errors along the beam center axis were extracted, and their FOM values ​​were calculated. The FOM values ​​of each dose along the beam center axis were compared between the two methods. Figure 2 As shown. Figure 2 The horizontal axis represents the distance from a certain location inside the human body in the voxel model to the skin surface, and the vertical axis represents the FOM value of the dose. The unmarked curves in the figure represent the FOM values ​​of boron dose, hydrogen dose, nitrogen dose, photon dose, and total dose calculated by the Ray Tracking method, while the marked curves represent the FOM values ​​of boron dose, hydrogen dose, nitrogen dose, photon dose, and total dose calculated by the method proposed in this invention. According to... Figure 2 On the beam center axis, the method proposed in this invention improves the calculation efficiency of boron dose, hydrogen dose, and photon dose to a certain extent. The calculation efficiency of hydrogen dose is significantly improved, while the calculation efficiency of nitrogen dose is slightly lower than that of Ray Tracking. However, the total dose FOM value is still significantly improved, which verifies the effectiveness of the method proposed in this invention.

[0070] In the description disclosed in this invention, it should be understood that the terms "center", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. They are only for the convenience of describing the invention and simplifying the description, and do not indicate or imply that the device or element referred to has a specific orientation, or is constructed and operated in a specific orientation. Therefore, they should not be construed as limiting the disclosure of this invention.

[0071] Unless the context otherwise requires, throughout the specification and claims, the term "comprising" is interpreted as open-ended and encompassing, meaning "including, but not limited to." In the description, terms such as "one embodiment," "some embodiments," "exemplary embodiment," "exemplary," or "some examples" are intended to indicate that a particular feature, structure, material, or characteristic associated with that embodiment or example is included in at least one embodiment or example disclosed in the invention. The illustrative representations of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics mentioned may be included in any suitable manner in any one or more embodiments or examples.

[0072] The terms "first" and "second" are used merely to distinguish different descriptive objects and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated; that is, they do not limit the position, order, priority, quantity, or content of the described objects. Therefore, a feature specified as "first" or "second" may explicitly or implicitly include one or more of that feature. In the description of the embodiments disclosed in this invention, unless otherwise stated, "a plurality of" means two or more.

[0073] In describing some embodiments, the term "connection" and its derivative expressions may be used. For example, the term "connection" may be used in describing some embodiments to indicate that two or more components have direct physical or electrical contact with each other. The embodiments disclosed herein are not necessarily limited to the content of this document.

[0074] "At least one of A, B, and C" has the same meaning as "at least one of A, B, or C," both including the following combinations of A, B, and C: only A, only B, only C, a combination of A and B, a combination of A and C, a combination of B and C, and a combination of A, B, and C. In the description of the embodiments of this application, unless otherwise stated, " / " means "or," for example, A / B can mean A or B; "A and / or B" includes the following three combinations: only A, only B, and a combination of A and B.

[0075] The use of “configured as” in this article implies an open and inclusive language that does not exclude the applicability to or configuration of devices to perform additional tasks or steps.

[0076] In addition, the use of “based on” implies openness and inclusivity, because processes, steps, calculations or other actions “based on” one or more of the stated conditions or values ​​may in practice be based on additional conditions or values ​​beyond those stated.

[0077] As used herein, “about” and “approximately” include the values ​​stated and the average values ​​within an acceptable range of deviation from a particular value, wherein the acceptable range of deviation is determined by a person skilled in the art taking into account the measurement under discussion and the error associated with the measurement of the particular quantity (i.e., the limitations of the measurement system).

[0078] This document describes exemplary embodiments with reference to cross-sectional views and / or plan views, which are idealized exemplary drawings. In the drawings, the thickness of layers and regions is enlarged for clarity. Therefore, variations in shape relative to the drawings are contemplated due to, for example, manufacturing techniques and / or tolerances. Thus, exemplary embodiments should not be construed as limited to the shapes of the regions shown herein, but rather include shape deviations due to, for example, manufacturing processes. For example, etched regions shown as rectangular would typically have curved features. Therefore, the regions shown in the drawings are schematic in nature, and their shapes are not intended to show the actual shapes of the regions of the device, nor are they intended to limit the scope of the exemplary embodiments.

[0079] The above description is merely a specific embodiment of the present invention, but the scope of protection disclosed in the present invention is not limited thereto. Any variations or substitutions conceived by those skilled in the art within the scope of the technology disclosed in the present invention should be included within the scope of protection disclosed in the present invention. Therefore, the scope of protection disclosed in the present invention should be determined by the scope of the claims.

Claims

1. A BNCT dose calculation method, characterized in that, Includes the following steps: S1. Construct a voxel model and initialize one or more particles in the simulation system, and obtain the initial position, energy and direction of motion of each particle respectively; S2. Based on the region where the particle is located and the energy of the particle, adaptively select the Monte Carlo particle tracking algorithm and flux counting method to count the dose; S3. Repeat step S2 until the particle leaves the simulation system or is killed by meeting the cutoff condition; S4. Cumulative calculation of dose distribution.

2. The BNCT dose calculation method according to claim 1, characterized in that, In step S1, the voxel model is a digital model, obtained by reading information from the medical image files of the case and converting it into a format that can be recognized by the Monte Carlo particle transport program.

3. The BNCT dose calculation method according to claim 1, characterized in that, In step S2, the implementation of the adaptive selection of the Monte Carlo particle tracking algorithm and the flux counting method includes: S21. Determine the region of each particle; S22. If the particle is located in the air region, Delta Tracking is used for particle tracking; if the particle is located inside the human body region, proceed to step S23. S23. If the energy of the particle is higher than the preset energy threshold E th If the particle tracking is positive, Ray Tracking is used; otherwise, Delta Tracking is used.

4. The BNCT dose calculation method according to claim 1, characterized in that, In step S3, the cutoff condition is any one of particle energy cutoff, particle weight cutoff, or time cutoff.

5. The BNCT dose calculation method according to claim 1, characterized in that, In step S4, the step of accumulating the dose distribution includes: calculating the photon equivalent dose of the voxel grid in the voxel model in each voxel grid, and accumulating it into the dose distribution of the voxel model.

6. The BNCT dose calculation method according to claim 3, characterized in that, The Delta Tracking includes: a cross-sectional value Σ based on a cross-sectional value greater than or equal to the maximum total macroscopic cross-section of all materials within the simulation system. maj The transport length of the particle is sampled, and after the particle collides, the macroscopic total cross section of the material corresponding to the voxel mesh where the particle is currently located is calculated based on Σ. maj The ratio is used to determine whether the collision is a real physical collision or a virtual collision in a probabilistic way.

7. The BNCT dose calculation method according to claim 3, characterized in that, The Ray Tracking includes: directly sampling the distance from the particle to the next physical collision point based on the macroscopic total cross section of the material corresponding to the voxel grid where the particle is currently located, calculating the distance from the particle to the next geometric boundary along the current direction, and taking the smaller of the two values ​​as the next track length of the particle.

8. An electronic device, characterized in that, Includes a processor and a memory, wherein the processor is coupled to the memory; The memory is used to store a computer program, and the processor is used to execute the computer program stored in the memory to implement the BNCT dose calculation method as described in any one of claims 1 to 7.

9. A computer-readable storage medium for storing a computer program, characterized in that, The computer program is executed by a processor to implement the BNCT dose calculation method as described in any one of claims 1 to 7.

Images