Neutron dose calculation method, device and equipment and readable storage medium

Abstract

The invention relates to a neutron dose calculation method, device and equipment and a readable storage medium, and belongs to the technical field of boron neutron capture therapy dose calculation.The method comprises the steps that an initial voxel grid model generated based on medical image information is obtained, the initial voxel grid model comprises multiple voxels, and the multiple voxels comprise air voxels; cutting all air voxels in the initial voxel grid model to generate a target voxel grid model; a plurality of source particles used for Monte Carlo simulation are obtained, the plurality of source particles are simulated in the target voxel grid model based on an event-based Monte Carlo algorithm, a neutron dose distribution result of the target voxel grid model is obtained, and the neutron dose distribution result comprises the neutron dose of each voxel in the target voxel grid model. According to the method, the air voxels are accurately cut, so that the efficiency and the accuracy of Monte Carlo simulation are improved, and the neutron dose can be calculated more quickly and accurately.

Description

Technical Field

[0001] The present invention relates to the technical field of boron neutron capture therapy dose calculation, and particularly to a method, device, equipment and readable storage medium for calculating neutron dose. Background Art

[0002] In boron neutron capture therapy (BNCT), in order to ensure that the dose at the tumor site of the patient reaches the prescribed dose and at the same time ensure that the dose in normal tissues does not exceed the limit, it is necessary to import medical images such as CT and MRI of the patient into the treatment planning system, and combine neutron source information for dose calculation and analysis.

[0003] BNCT dose calculation is a photon-coupled transport problem in complex geometries, and the generated dose mainly includes 10 B(n,α) 6 boron dose generated by the Li reaction, 1 H(n,n’) 1 hydrogen dose generated by the H reaction, 14 N(n,p) 14 nitrogen dose generated by the C reaction, photon dose generated by photons in the beam and secondary photons.

[0004] To ensure the accuracy of dose calculation, existing BNCT dose calculations usually use Monte Carlo programs as the core of dose calculation in the treatment planning system. The geometric model of the Monte Carlo program is directly generated from the patient's CT and MRI, generally including 20 million to 50 million voxel grids. To ensure the convergence of dose statistics results in the voxel grids, a large number of particles need to be simulated, which makes the dose calculation process extremely time-consuming. Although the form of multi-core parallelism can accelerate the calculation, due to the complexity of the problem itself and the huge amount of calculation, each dose calculation still takes several hours, and the calculation speed is still too slow, seriously restricting the formulation efficiency of BNCT treatment plans. In the simulation process of a single source particle, a large number of judgment statements are also involved. Directly applying the algorithm to GPU calculation will lead to serious thread divergence, and the calculation speed is even slower than that of the CPU platform. Summary of the Invention

[0005] The technical problem to be solved by the present invention is to provide a method, device, equipment and readable storage medium for calculating neutron dose, aiming to solve at least one of the above technical problems.

[0006] The technical solution of the present invention to solve the above technical problems is as follows:

[0007] In the first aspect, the present application provides a method, device, equipment and readable storage medium for calculating neutron dose, and adopts the following technical solutions:

[0008] A method for calculating neutron dose includes:

[0009] Obtain an initial voxel grid model generated based on medical image information. The initial voxel grid model includes a plurality of voxels, and the plurality of voxels includes air voxels, and the air voxels represent the volume elements occupied by air;

[0010] Crop all the air voxels in the initial voxel grid model to generate a target voxel grid model;

[0011] Obtain a plurality of source particles for Monte Carlo simulation, and based on the event-based Monte Carlo algorithm, simulate the plurality of source particles in the target voxel grid model to obtain the neutron dose distribution result of the target voxel grid model. The neutron dose distribution result includes the neutron dose of each voxel in the target voxel grid model.

[0012] The beneficial effects of the present invention are as follows: The method first crops all the air voxels in the initial voxel grid model to generate a target voxel grid model, effectively reducing the simulation of irrelevant regions and reducing the occupancy of GPU video memory. Then, based on the event-based Monte Carlo algorithm, the same calculation processes of different source particles are packed, and the plurality of source particles are parallelly simulated in the target voxel grid model, avoiding the thread divergence phenomenon caused by a large number of judgment statements and conditional branches when the traditional history-based algorithm is applied on the GPU because each thread independently simulates the complete history of a source particle, which seriously affects the parallel computing efficiency of the GPU. The computing efficiency and parallel performance are greatly improved, thereby accelerating the speed of BNCT dose calculation and improving the accuracy and reliability of dose calculation.

[0013] On the basis of the above technical solution, the present invention can also be improved as follows.

[0014] Further, the step of cropping all the air voxels in the initial voxel grid model to generate a target voxel grid model includes:

[0015] Establish a coordinate system of the initial voxel grid model, obtain a scan result for the initial voxel grid model, and based on the scan result, determine the boundary points of non-air voxels on the X-axis of the coordinate system and the boundary points of non-air voxels on the Y-axis in each layer of voxel grid;

[0016] Based on the boundary points of non-air voxels on the X-axis and the boundary points of non-air voxels on the Y-axis in all layers of voxel grids, determine boundary information;

[0017] Construct a bounding box according to the boundary information;

[0018] Based on the bounding box, all air voxels in the initial voxel grid model are cropped to generate a target voxel grid model.

[0019] The beneficial effect of adopting the above further solution is that all air voxels in the initial voxel grid model are cropped to generate a target voxel grid model. Specifically, first, the initial voxel grid model is scanned, and based on the scanning results, the boundary points of non-air voxels on the X-axis and the boundary points of non-air voxels on the Y-axis in each layer of voxel grid are determined. Then, based on the boundary points of non-air voxels on the X-axis in all voxel grids and the boundary points of non-air voxels on the Y-axis in all voxel grids, the boundary information is determined. Next, according to the boundary information, a bounding box is constructed. Finally, based on the bounding box, all air voxels in the initial voxel grid model are cropped to generate a target voxel grid model. This process effectively reduces irrelevant air voxels, reduces the GPU video memory occupancy of the geometric model, dose matrix, and error matrix, thereby improving the program running speed and enhancing the efficiency of BNCT dose calculation.

[0020] Further, determining the boundary information based on the boundary points of non-air voxels on the X-axis and the boundary points of non-air voxels on the Y-axis in the coordinate system in all layers of voxel grids includes:

[0021] Merge the boundary points of non-air voxels on the X-axis in all layers of voxel grids in the coordinate system to obtain the maximum boundary point and the minimum boundary point of non-air voxels in the initial voxel grid model in the X-axis direction of the coordinate system;

[0022] Merge the boundary points of non-air voxels on the Y-axis in all layers of voxel grids in the coordinate system to obtain the maximum boundary point and the minimum boundary point of non-air voxels in the initial voxel grid model in the Y-axis direction of the coordinate system;

[0023] Based on the maximum boundary point and the minimum boundary point of non-air voxels in the initial voxel grid model in the X-axis direction and the Y-axis direction of the coordinate system respectively, determine the boundary information.

[0024] The beneficial effect of adopting the above further solution is: The boundary information of non-air voxels in the initial voxel grid model can be accurately determined. Specifically, by merging the boundary points of non-air voxels on the X-axis and the Y-axis in each layer of voxel grids, the maximum and minimum boundary points of non-air voxels in the entire model are obtained, thereby effectively reducing unnecessary air voxels, reducing the GPU video memory occupancy required for subsequent dose calculations, and improving the calculation efficiency.

[0025] Further, in the event-based Monte Carlo algorithm, simulating the multiple source particles in the target voxel grid model to obtain the neutron dose distribution result of the target voxel grid model includes:

[0026] Step S31: Select a set number of source particles from the multiple source particles;

[0027] Step S32: Generate a batch of source particles to be simulated based on the selected source particles;

[0028] Step S33: Based on multiple threads of the GPU itself, perform parallel simulation on all source particles in the batch of source particles to be simulated in the target voxel grid model to obtain the local neutron dose corresponding to each thread;

[0029] Step S34: Determine whether there are unselected source particles;

[0030] Step S35: If there are unselected source particles, select a set number of source particles from the unselected multiple source particles, and execute Steps S32 to S35 until all source particles are simulated in the target voxel grid model. Based on the multiple local neutron doses corresponding to each thread, determine the neutron dose distribution result of the target voxel grid model.

[0031] The beneficial effect of adopting the above further solution is as follows: Based on the event-based Monte Carlo algorithm, simulate multiple source particles in the target voxel grid model. The specific steps include: selecting a set number of source particles and generating a batch of source particles to be simulated; using multiple threads of the GPU itself to perform parallel simulation on these source particles to obtain the local neutron dose corresponding to each thread; determining whether there are still unselected source particles. If so, continue to select new source particles and repeat the above process until all source particles are simulated. Finally, based on the local neutron doses corresponding to each thread, comprehensively obtain the overall neutron dose distribution result of the target voxel grid model. This solution effectively alleviates the GPU thread divergence phenomenon, improves the GPU utilization rate, and significantly speeds up the BNCT dose calculation speed.

[0032] Further, based on multiple threads of the GPU itself, performing parallel simulation on all source particles in the batch of source particles to be simulated in the target voxel grid model to obtain the local neutron dose corresponding to each thread includes:

[0033] Based on a preset particle allocation rule, allocate each source particle in the batch of source particles to be simulated to the threads of the GPU itself, where each thread is allocated a selected particle;

[0034] Based on a preset simulation process, all source particles in the batch of source particles to be simulated are subjected to parallel simulation in the target voxel grid model until the particles simulated by each thread are absorbed by the target voxel model, and the energy deposition of each particle in each voxel is calculated to obtain the local neutron dose corresponding to each thread.

[0035] The beneficial effects of adopting the above further solution are as follows: It realizes the parallel simulation of all source particles in the batch of source particles to be simulated in the target voxel grid model, calculates the energy deposition of each source particle in each voxel, and thus obtains the local neutron dose corresponding to each thread. This method makes full use of the parallel processing ability of GPU multi-threads, effectively improves the speed and efficiency of Monte Carlo simulation, reduces the calculation time, and is applicable to neutron dose calculation under a large-scale voxel grid.

[0036] The parallel simulation of all source particles in the batch of source particles to be simulated in the target voxel grid model based on the preset simulation process includes:

[0037] Step S331: For each thread, calculate the reaction cross-section of the source particles corresponding to this thread in the target voxel model, and the value of the reaction cross-section characterizes the probability of nuclear reaction between the source particles and the target nuclei in the target voxel model;

[0038] Step S332: For each thread, based on the transport length sampling kernel and the reaction cross-section of the source particles corresponding to this thread, perform transport length sampling on the source particles corresponding to this thread in the target voxel model to obtain the transport length of the source particles, and the value of the transport length characterizes the distance that the source particles move in the target voxel model;

[0039] Step S333: For each thread, calculate the nearest boundary distance of the source particles corresponding to this thread along the movement direction to the voxel where they are located based on the boundary distance calculation kernel;

[0040] Step S334: For each thread, perform absorption processing on the source particles corresponding to this thread in the target voxel model based on the absorption reaction kernel to obtain the absorption state of the source particles corresponding to this thread;

[0041] Step S335: For each thread, perform scattering processing on the source particles corresponding to this thread in the target voxel model based on the scattering reaction kernel to obtain the scattering information of the source particles corresponding to this thread;

[0042] Step S336: For each of the said threads, based on the KERMA calculation kernel, calculate the KERMA values of multiple elements in the voxel where the source particles corresponding to the thread are located. The KERMA value characterizes the initial kinetic energy of charged particles ionized per unit mass of the target material by the source particles. The multiple elements include boron, hydrogen, and nitrogen elements.

[0043] Step S337: For each of the said threads, based on the dose counting kernel and the KERMA value of the source particles corresponding to the thread, perform dose calculation on the source particles corresponding to the thread to obtain the local neutron dose corresponding to the thread.

[0044] Step S338: For each of the said threads, based on the transport length, nearest boundary distance, absorption state, and scattering information of the source particles corresponding to the thread, determine whether there are source particles in the target voxel model.

[0045] Step S339: For each of the said threads, if there are still source particles in the target voxel model, execute Steps S332 to S339 until there are no particles in the target voxel model, and the simulation process ends.

[0046] The beneficial effects of adopting the above further solution are as follows: By using the event-based Monte Carlo algorithm, parallel simulation of multiple source particles in the target voxel grid model is carried out, which greatly alleviates the thread divergence phenomenon of the GPU and improves the utilization rate of the GPU. Specifically, the same calculation processes of different source particles are packaged, and the same processes such as transport length sampling, implicit capture, and energy-angle distribution sampling of particles in different states are merged, and are respectively written into different kernel functions. By calling each kernel function, the reaction cross-section, transport length sampling, nearest boundary distance calculation, absorption processing, scattering processing, KERMA value calculation, and dose counting of source particles are calculated for each thread, realizing the efficient parallel simulation of the same transport process of different particles by multiple threads. This not only speeds up the dose calculation speed but also ensures the convergence of the dose statistical results in the voxel grid, thus significantly improving the overall efficiency of BNCT dose calculation.

[0047] In the second aspect, the present application provides an application of a method for calculating neutron dose, adopting the following technical solution:

[0048] An application of a method for calculating neutron dose in boron neutron capture therapy.

[0049] In the third aspect, the present application provides a device for calculating neutron dose, adopting the following technical solution:

[0050] A device for calculating neutron dose, comprising:

[0051] An acquisition module, configured to acquire an initial voxel grid model generated based on medical image information, where the initial voxel grid model includes a plurality of voxels, and the plurality of voxels include air voxels, and the air voxels represent volume elements occupied by air;

[0052] A cutting module, configured to cut all air voxels in the initial voxel grid model to generate a target voxel grid model;

[0053] A simulation module, configured to acquire a plurality of source particles for Monte Carlo simulation, and based on the event-based Monte Carlo algorithm, simulate the plurality of source particles in the target voxel grid model to obtain a neutron dose distribution result of the target voxel grid model, where the neutron dose distribution result includes the neutron dose of each voxel in the target voxel grid model.

[0054] In a fourth aspect, the present application provides an electronic device, adopting the following technical solution:

[0055] An electronic device includes a memory and a processor, and a computer program capable of being loaded and executed by the processor for the calculation method of neutron dose according to any one of the first aspects is stored on the memory.

[0056] In a fifth aspect, the present application provides a computer-readable storage medium, adopting the following technical solution:

[0057] A computer-readable storage medium stores a computer program capable of being loaded and executed by the processor for the calculation method of neutron dose according to any one of the first aspects.

[0058] Additional aspects and advantages of the present application will be given in part in the following description, and these will become obvious from the following description, or can be understood through the practice of the present application. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1 It is a schematic flowchart of a calculation method of neutron dose provided by an embodiment of the present invention;

[0060] Figure 2 It is a schematic flowchart of a method for generating a target voxel grid model provided by an embodiment of the present invention;

[0061] Figure 3 It is a schematic structural diagram of an event-based Monte Carlo algorithm provided by an embodiment of the present invention;

[0062] Figure 4 It is a schematic structural diagram of another event-based Monte Carlo algorithm provided by an embodiment of the present invention;

[0063] Figure 5A schematic diagram of the structure of a neutron dose calculation device provided by one embodiment of the present invention;

[0064] Figure 6 A schematic diagram of the structure of an electronic device provided by an embodiment of the present invention. DETAILED DESCRIPTION

[0065] To make the purpose, technical solutions and advantages of the embodiments of the present application clearer, the following is a Figures 1 to 6 , the technical solutions in the embodiments of the present application are clearly and completely described. Obviously, the described embodiments are part of the embodiments of the present application, not all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present application.

[0066] In addition, the term "and / or" in this article is only a description of the association relationship of associated objects, indicating that there can be three relationships. For example, A and / or B can represent: A exists alone, A and B exist at the same time, and B exists alone. In addition, the character " / " in this article, unless otherwise specified, generally means that the associated objects before and after are in an "or" relationship.

[0067] like Figure 1 As shown, an embodiment of the present application provides a method for calculating a neutron dose, and the application of the method for calculating a neutron dose in boron neutron capture therapy. A method for calculating a neutron dose mainly includes steps S1 to S3:

[0068] Step S1, obtaining an initial voxel grid model generated based on medical image information, wherein the initial voxel grid model includes a plurality of voxels, wherein the plurality of voxels include air voxels, and the air voxels represent volume elements occupied by air;

[0069] In the embodiment of the present application, first, key information such as tumor location, normal tissue boundary, etc. is extracted from medical images such as CT or MRI, and then the extracted key information is converted into a three-dimensional voxel grid model, where each voxel represents a small cubic space and contains material properties and density information. The OpenCV library or the VTK library can be used for image processing and data conversion.

[0070] Medical imaging can also support other types of data such as PET (positron emission tomography) and SPECT (single photon emission computed tomography). These imaging data provide richer physiological and metabolic information, which helps to locate tumors more accurately and evaluate the effectiveness of treatment.

[0071] Step S2, cutting all air voxels in the initial voxel grid model to generate a target voxel grid model;

[0072] In the embodiments of the present application, as Figure 2 shown, specifically, step S2 includes the following sub-steps:

[0073] Step S21, establish a coordinate system for the initial voxel grid model, obtain a scanning result for the initial voxel grid model, and based on the scanning result, determine the boundary points of non-air voxels on the X-axis of the coordinate system and the boundary points of non-air voxels on the Y-axis in each layer of the voxel grid;

[0074] In the embodiments of the present application, during the scanning process, the device or software will traverse each voxel in the model and record its attributes such as position and density. In the scanning result, air voxels and non-air voxels are distinguished according to the density attribute of the voxels. Generally, the density of air voxels is relatively low. When the density of any voxel is lower than a preset density threshold, the voxel is determined to be an air voxel.

[0075] In the embodiments of the present application, starting from the top layer of the initial voxel grid model, traverse the voxel grid layer by layer downward. For each non-air voxel in the current layer, check its adjacent voxels in the X-axis direction. If there is an air voxel among the adjacent voxels and there is no non-air voxel on the other side of the current voxel in the X-axis direction, then the current voxel is the boundary point of non-air voxels on the X-axis. Similarly, the boundary points of non-air voxels on the Y-axis are also determined in the above manner and will not be elaborated here.

[0076] Step S22, based on the boundary points of non-air voxels on the X-axis of the coordinate system and the boundary points of non-air voxels on the Y-axis in the voxel grids of all layers, determine the boundary information;

[0077] In the embodiments of the present application, step S22 includes:

[0078] Merge the boundary points of non-air voxels on the X-axis of the coordinate system in the voxel grids of all layers to obtain the maximum boundary point and the minimum boundary point of non-air voxels in the initial voxel grid model in the X-axis direction of the coordinate system;

[0079] Merge the boundary points of non-air voxels on the Y-axis of the coordinate system in the voxel grids of all layers to obtain the maximum boundary point and the minimum boundary point of non-air voxels in the initial voxel grid model in the Y-axis direction of the coordinate system;

[0080] Based on the maximum boundary point and the minimum boundary point of non-air voxels in the initial voxel grid model in the X-axis direction and the Y-axis direction of the coordinate system respectively, determine the boundary information.

[0081] In the embodiments of the present application, the boundary information describes the spatial range of air voxels in the initial voxel grid model in the X-axis and Y-axis directions.

[0082] Step S23: Construct a bounding box according to the boundary information.

[0083] In the embodiments of the present application, an appropriate software tool or programming language (such as MATLAB, Python, etc.) is used to construct a bounding box according to the determined dimensions.

[0084] Step S24: Based on the bounding box, cut all the air voxels in the initial voxel grid model to generate a target voxel grid model.

[0085] In the embodiments of the present application, for each voxel, check whether its position is inside the bounding box, and cut the voxels outside the bounding box to generate a target voxel grid model.

[0086] Step S3: Obtain a plurality of source particles for Monte Carlo simulation, and based on the event-based Monte Carlo algorithm, simulate the plurality of source particles in the target voxel grid model to obtain the neutron dose distribution result of the target voxel grid model. The neutron dose distribution result includes the neutron dose of each voxel in the target voxel grid model, and the neutron dose includes boron dose, hydrogen dose, and nitrogen dose.

[0087] In the embodiments of the present application, as Figure 3 and Figure 4 shown, the process of simulating the plurality of source particles in the target voxel grid model based on the event-based Monte Carlo algorithm to obtain the neutron dose distribution result of the target voxel grid model includes:

[0088] Step S31: Select a set number of source particles from the plurality of source particles;

[0089] In the embodiments of the present application, the set number is determined according to the multiple threads of the GPU.

[0090] Step S32: Generate a batch of source particles to be simulated based on the selected source particles;

[0091] Step S33: Based on the multiple threads of the GPU itself, perform parallel simulation on all the source particles in the batch of source particles to be simulated in the target voxel grid model to obtain the local neutron dose corresponding to each thread;

[0092] In the embodiments of the present application, all the source particles in the batch of source particles to be simulated are assigned to the multiple threads of the GPU for parallel simulation. Each thread simulates the movement trajectory and interaction process of the source particles assigned to it in the target voxel grid model, and obtains the local neutron dose corresponding to each thread by counting the number of neutrons and neutron energy in each voxel.

[0093] Step S34, determine whether there are unselected source particles, and check whether there are still source particles that have not been simulated for the next round of simulation;

[0094] Step S35, if there are unselected source particles, select a set number of source particles from the multiple unselected source particles, and execute Steps S32 to S35 until all source particles have been simulated in the target voxel grid model, ensuring that all source particles have been simulated in the target voxel grid model, and determining the neutron dose distribution result of the target voxel grid model based on the multiple local neutron doses corresponding to each thread, so as to obtain an accurate neutron dose distribution result.

[0095] In the embodiment of the present application, the local neutron doses obtained in each round of simulation are accumulated to obtain the final neutron dose distribution result.

[0096] In the above embodiment, based on multiple threads of the GPU itself, all source particles in the batch of source particles to be simulated are parallelly simulated in the target voxel grid model to obtain the local neutron dose corresponding to each thread, including:

[0097] Based on a preset particle distribution rule, each source particle in the batch of source particles to be simulated is respectively assigned to a thread of the GPU itself, where each thread is assigned a selected source particle;

[0098] Based on a preset simulation process, all source particles in the batch of source particles to be simulated are parallelly simulated in the target voxel grid model, and the energy deposition of each source particle in each voxel is calculated to obtain the local neutron dose corresponding to each thread.

[0099] Specifically, based on the preset simulation process, parallelly simulating all source particles in the batch of source particles to be simulated in the target voxel grid model includes:

[0100] Step S331, for each thread, calculate the reaction cross-section of the source particle corresponding to the thread in the target voxel model, and the value of the reaction cross-section characterizes the probability of the source particle undergoing a nuclear reaction with the target nucleus in the target voxel model;

[0101] Step S332, for each thread, based on the transport length sampling kernel and the reaction cross-section of the source particle corresponding to the thread, perform transport length sampling on the source particle corresponding to the thread in the target voxel model to obtain the transport length of the source particle, and the value of the transport length characterizes the distance that the source particle moves in the target voxel model;

[0102] Step S333: For each of the threads, calculate the kernel based on the boundary distance, and calculate the nearest boundary distance of the source particle corresponding to the thread in the target voxel model;

[0103] Step S334: For each of the threads, based on the absorption reaction kernel, perform absorption processing on the source particle corresponding to the thread in the target voxel model to obtain the absorption state of the source particle corresponding to the thread;

[0104] Step S335: For each of the threads, based on the scattering reaction kernel, perform scattering processing on the source particle corresponding to the thread in the target voxel model to obtain the scattering information of the source particle corresponding to the thread;

[0105] Step S336: For each of the threads, based on the KERMA calculation kernel, calculate the KERMA values of multiple elements in the voxel where the source particle corresponding to the thread is located. The KERMA value characterizes the initial kinetic energy of charged particles ionized per unit mass of the target material by the source particle. The multiple elements include boron, hydrogen, and nitrogen;

[0106] Step S337: For each of the threads, based on the dose counting kernel and the KERMA value of the source particle corresponding to the thread, perform dose calculation on the source particle corresponding to the thread to obtain the local neutron dose corresponding to the thread;

[0107] Step S338: For each of the threads, based on the transport length, nearest boundary distance, absorption state, and scattering information of the source particle corresponding to the thread, determine whether there is a source particle in the target voxel model;

[0108] Step S339: For each of the threads, if there is a particle in the target voxel model, execute Steps S332 to S339 until there is no particle in the target voxel model, and the simulation process ends.

[0109] The same calculation processes for different source particles are packaged, and the same processes such as sampling the transport length, implicit capture, and sampling the energy-angle distribution of particles in different states are merged and written into different kernel functions respectively. By calling each kernel function, the reaction cross-section, transport length sampling, nearest boundary distance calculation, absorption processing, scattering processing, KERMA value calculation, and dose counting of source particles are calculated in parallel for each thread, realizing the efficient parallel simulation of the same transport process of different particles by multiple threads.

[0110] Modify the simulation of multiple threads for the complete independent histories of multiple particles to the simulation of multiple threads for the same calculation process of multiple particles, which can greatly alleviate the thread divergence phenomenon of the GPU, improve the utilization rate of the GPU, and accelerate the dose calculation of BNCT. It not only speeds up the dose calculation, but also ensures the convergence of the dose statistical results in the voxel grid, thus significantly improving the overall efficiency of BNCT dose calculation.

[0111] This method first cuts all the air voxels in the initial voxel grid model to generate a target voxel grid model, effectively reducing the simulation of irrelevant regions and reducing the occupancy of GPU video memory. Then, based on the event-based Monte Carlo algorithm, the same calculation processes of different source particles are packed, and multiple source particles are simulated in parallel in the target voxel grid model, avoiding the thread divergence phenomenon caused by a large number of judgment statements and conditional branches when the traditional history-based algorithm is applied on the GPU, where each thread independently simulates the complete history of a source particle, which seriously affects the parallel computing efficiency of the GPU. It greatly improves the computing efficiency and parallel performance, thus speeding up the BNCT dose calculation and improving the accuracy and reliability of the dose calculation.

[0112] Figure 5 The structural schematic diagram of a neutron dose calculation device 200 is shown.

[0113] As Figure 5 shown, a neutron dose calculation device 200 mainly includes:

[0114] An acquisition module 201, configured to acquire an initial voxel grid model generated based on medical image information, where the initial voxel grid model includes multiple voxels, and the multiple voxels include air voxels, and the air voxels represent the volume elements occupied by air;

[0115] A cutting module 202, configured to cut all the air voxels in the initial voxel grid model to generate a target voxel grid model;

[0116] A simulation module 203, configured to acquire multiple source particles for Monte Carlo simulation, and based on the event-based Monte Carlo algorithm, simulate the multiple source particles in the target voxel grid model to obtain the neutron dose distribution result of the target voxel grid model, where the neutron dose distribution result includes the neutron dose of each voxel in the target voxel grid model, and the neutron dose includes boron dose, hydrogen dose, and nitrogen dose.

[0117] Optionally, the cutting module 202 includes:

[0118] A scanning sub-module that establishes the coordinate system of the initial voxel grid model, obtains the scanning result for the initial voxel grid model, and based on the scanning result, determines the boundary points of non-air voxels on the X-axis of the coordinate system and the boundary points of non-air voxels on the Y-axis in each layer of voxel grids; A determining sub-module that is used to determine boundary information based on the boundary points of non-air voxels on the X-axis and the boundary points of non-air voxels on the Y-axis in all voxel grids;

[0119] A constructing sub-module that is used to construct a bounding box according to the boundary information;

[0120] A cutting sub-module that is used to cut all air voxels in the initial voxel grid model based on the bounding box to generate a target voxel grid model.

[0121] Optionally, the determining sub-module is specifically used for:

[0122] Merge the boundary points of non-air voxels on the X-axis of the coordinate system in the voxel grids of all layers to obtain the maximum boundary point and the minimum boundary point of non-air voxels in the initial voxel grid model in the X-axis direction of the coordinate system;

[0123] Merge the boundary points of non-air voxels on the Y-axis of the coordinate system in the voxel grids of all layers to obtain the maximum boundary point and the minimum boundary point of non-air voxels in the initial voxel grid model in the Y-axis direction of the coordinate system;

[0124] Determine the boundary information based on the maximum boundary point and the minimum boundary point of non-air voxels in the initial voxel grid model in the X-axis direction and the Y-axis direction of the coordinate system.

[0125] Optionally, the simulation module 203 is specifically used to perform the following steps:

[0126] Step S31, select a set number of source particles from the multiple source particles;

[0127] Step S32, generate a batch of source particles to be simulated based on the selected source particles;

[0128] Step S33, based on multiple threads of the GPU itself, perform parallel simulation of all source particles in the batch of source particles to be simulated in the target voxel grid model to obtain the local neutron dose corresponding to each thread;

[0129] Step S34, determine whether there are unselected source particles;

[0130] Step S35. If there are unselected source particles, select a set number of source particles from the multiple unselected source particles, and execute Steps S32 to S35 until all source particles are simulated in the target voxel grid model. Based on the multiple local neutron doses corresponding to each thread, determine the neutron dose distribution result of the target voxel grid model.

[0131] Optionally, based on multiple threads of the GPU itself, perform parallel simulation of all source particles in the batch of source particles to be simulated in the target voxel grid model to obtain the local neutron dose corresponding to each thread, including:

[0132] Based on a preset particle allocation rule, allocate each source particle in the batch of source particles to be simulated to the threads of the GPU itself, where each thread is allocated a selected source particle.

[0133] Based on a preset simulation process, perform parallel simulation of all source particles in the batch of source particles to be simulated in the target voxel grid model, and calculate the energy deposition of each source particle in each voxel to obtain the local neutron dose corresponding to each thread.

[0134] Optionally, based on a preset simulation process, performing parallel simulation of all source particles in the batch of source particles to be simulated in the target voxel grid model includes:

[0135] Step S331. For each thread, calculate the reaction cross-section of the source particle corresponding to the thread in the target voxel model, and the value of the reaction cross-section characterizes the probability of the source particle undergoing a nuclear reaction with the target nucleus in the target voxel model.

[0136] Step S332. For each thread, based on the transport length sampling kernel and the reaction cross-section of the source particle corresponding to the thread, perform transport length sampling of the source particle corresponding to the thread in the target voxel model to obtain the transport length of the source particle, and the value of the transport length characterizes the distance that the source particle moves in the target voxel model.

[0137] Step S333. For each thread, calculate the nearest boundary distance of the source particle corresponding to the thread in the target voxel model based on the boundary distance calculation kernel.

[0138] Step S334. For each thread, perform absorption processing on the source particle corresponding to the thread in the target voxel model based on the absorption reaction kernel to obtain the absorption state of the source particle corresponding to the thread.

[0139] Step S335: For each of the threads, based on the scattering reaction kernel, perform scattering processing on the source particles corresponding to the thread in the target voxel model to obtain the scattering information of the source particles corresponding to the thread;

[0140] Step S336: For each of the threads, based on the KERMA calculation kernel, calculate the KERMA values of multiple elements in the voxel where the source particles corresponding to the thread are located. The KERMA value characterizes the initial kinetic energy of the charged particles ionized per unit mass of the target material by the source particles;

[0141] Step S337: For each of the threads, based on the dose counting kernel and the KERMA value of the source particles corresponding to the thread, perform dose calculation on the source particles corresponding to the thread to obtain the local neutron dose corresponding to the thread;

[0142] Step S338: For each of the threads, based on the transport length, the distance to the nearest boundary, the absorption state, and the scattering information of the source particles corresponding to the thread, determine whether there are source particles in the target voxel model;

[0143] Step S339: For each of the threads, if there are still source particles in the target voxel model, execute steps S332 to S339 until there are no particles in the target voxel model, and the simulation process ends.

[0144] In one example, the modules in any of the above devices can be one or more integrated circuits configured to implement the above methods. For example: one or more application specific integrated circuits (ASICs), or, one or more digital signal processors (DSPs), or, one or more field programmable gate arrays (FPGAs), or a combination of at least two of these integrated circuit forms.

[0145] Again, when the modules in the device can be implemented in the form of a processing element scheduler, the processing element can be a general-purpose processor, such as a central processing unit (CPU) or other processors that can call programs. Again, these modules can be integrated together to be implemented in the form of a system-on-a-chip (SOC).

[0146] In this application, names may be assigned to various objects such as various messages / information / devices / network elements / systems / devices / actions / operations / processes / concepts, etc. It can be understood that these specific names do not constitute a limitation on the relevant objects, and the assigned names may change with factors such as the scenario, context, or usage habits. The understanding of the technical meaning of the technical terms in this application should be mainly determined from the functions and technical effects embodied / executed in the technical solution.

[0147] Those skilled in the art can clearly understand that for the convenience and simplicity of description, the specific working processes of the systems, devices, and modules described above can refer to the corresponding processes in the foregoing method embodiments, and will not be elaborated herein.

[0148] Those of ordinary skill in the art can realize that the modules and algorithm steps of each example described in combination with the embodiments disclosed herein can be implemented by electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are executed in a hardware or software manner depends on the specific application and design constraints of the technical solution. Professional technicians can use different methods to implement the described functions for each specific application, but such implementation should not be considered to exceed the scope of this application.

[0149] Figure 6 It is a block diagram of the structure of an electronic device 300 according to an embodiment of this application.

[0150] As Figure 6 shown, the electronic device 300 includes a processor 301 and a memory 302, and may further include one or more of an information input / output (I / O) interface 303, a communication component 304, and a communication bus 305.

[0151] Among them, the processor 301 is used to control the overall operation of the electronic device 300 to complete all or part of the steps in the above-mentioned neutron dose calculation method; the memory 302 is used to store various types of data to support the operation of the electronic device 300. These data may include, for example, instructions for any application or method operating on the electronic device 300, as well as application-related data. The memory 302 can be implemented by any type of volatile or non-volatile storage device or a combination thereof, such as Static Random Access Memory (SRAM), Electrically Erasable Programmable Read-Only Memory (EEPROM), Erasable Programmable Read-Only Memory (EPROM), Programmable Read-Only Memory (PROM), Read-Only Memory (ROM), magnetic memory, flash memory, one or more of a magnetic disk or an optical disc.

[0152] The I / O interface 303 provides an interface between the processor 301 and other interface modules. The above-mentioned other interface modules can be a keyboard, a mouse, buttons, etc. These buttons can be virtual buttons or physical buttons. The communication component 304 is used to test the wired or wireless communication between the electronic device 300 and other devices. Wireless communication, such as Wi-Fi, Bluetooth, Near Field Communication (NFC), 2G, 3G or 4G, or a combination of one or more of them. Therefore, the corresponding communication component 304 can include: a Wi-Fi component, a Bluetooth component, an NFC component.

[0153] The communication bus 305 may include a path for transmitting information between the above components. The communication bus 305 can be a PCI (Peripheral Component Interconnect) bus or an EISA (Extended Industry Standard Architecture) bus, etc. The communication bus 305 can be divided into an address bus, a data bus, a control bus, etc.

[0154] The electronic device 300 can be implemented by one or more application specific integrated circuits (ASICs), digital signal processors (DSPs), digital signal processing devices (DSPDs), programmable logic devices (PLDs), field programmable gate arrays (FPGAs), controllers, microcontrollers, microprocessors or other electronic components, and is used to execute the neutron dose calculation method given in the above embodiments.

[0155] The computer-readable storage medium provided by the embodiments of the present application will be introduced below. The computer-readable storage medium described below can be correspondingly referred to the neutron dose calculation method described above.

[0156] The present application also provides a computer-readable storage medium, on which a computer program is stored. When the computer program is executed by a processor, the steps of the above neutron dose calculation method are implemented.

[0157] The computer-readable storage medium may include: various media that can store program codes, such as USB flash drives, mobile hard disks, read-only memories (ROMs), random access memories (RAMs), magnetic disks or optical discs.

[0158] The term "including", "comprising" or any other variant thereof is intended to cover non-exclusive inclusion, such that a process, method, article or device including a series of elements includes not only those elements but also other elements not expressly listed, or also includes elements inherent to such process, method, article or device.

[0159] The above description is only the preferred embodiments of the present application and the description of the applied technical principles. Those skilled in the art should understand that the scope of the application involved in the present application is not limited to the technical solutions formed by the specific combination of the above technical features, and should also cover other technical solutions formed by any combination of the above technical features or their equivalent features without departing from the foregoing application concept. For example, the technical solutions formed by mutually replacing the above features with the technical features (but not limited to) having similar functions applied in the present application.

Claims

1. A method for calculating neutron dose, characterized in that, Including: Obtain an initial voxel grid model generated based on medical image information, where the initial voxel grid model includes a plurality of voxels, and the plurality of voxels include air voxels, and the air voxels represent the volume elements occupied by air; Crop all the air voxels in the initial voxel grid model to generate a target voxel grid model; Obtain a plurality of source particles for Monte Carlo simulation, and based on the event-based Monte Carlo algorithm, simulate the plurality of source particles in the target voxel grid model to obtain the neutron dose distribution result of the target voxel grid model, where the neutron dose distribution result includes the neutron dose of each voxel in the target voxel grid model.

2. The calculation method of neutron dose according to claim 1, characterized in that, The step of cropping all the air voxels in the initial voxel grid model to generate a target voxel grid model includes: Establish a coordinate system for the initial voxel grid model, obtain a scan result for the initial voxel grid model, and based on the scan result, determine the boundary points of non-air voxels on the X-axis of the coordinate system and the boundary points of non-air voxels on the Y-axis in each layer of the voxel grid; Based on the boundary points of non-air voxels on the X-axis of the coordinate system and the boundary points of non-air voxels on the Y-axis in all layers of the voxel grid, determine the boundary information; Construct a bounding box according to the boundary information; Based on the bounding box, crop all the air voxels in the initial voxel grid model to generate a target voxel grid model.

3. A method for calculating neutron dose according to claim 2, characterized in that, The step of determining the boundary information based on the boundary points of non-air voxels on the X-axis of the coordinate system and the boundary points of non-air voxels on the Y-axis in all layers of the voxel grid includes: Merge the boundary points of non-air voxels on the X-axis of the coordinate system in all layers of the voxel grid to obtain the maximum boundary point and the minimum boundary point of non-air voxels in the initial voxel grid model in the X-axis direction of the coordinate system; Merge the boundary points of non-air voxels on the Y-axis of the coordinate system in all layers of the voxel grid to obtain the maximum boundary point and the minimum boundary point of non-air voxels in the initial voxel grid model in the Y-axis direction of the coordinate system; Based on the maximum boundary point and the minimum boundary point of non-air voxels in the initial voxel grid model in the X-axis direction and the Y-axis direction of the coordinate system respectively, determine the boundary information.

4. A method for calculating neutron dose according to claim 1, characterized in that, The step of simulating the plurality of source particles in the target voxel grid model based on the event-based Monte Carlo algorithm to obtain the neutron dose distribution result of the target voxel grid model includes: Step S31, select a set number of source particles from the plurality of source particles; Step S32, generate a batch of source particles to be simulated based on the selected source particles; Step S33, based on multiple threads of the GPU itself, perform parallel simulation of all the source particles in the batch of source particles to be simulated in the target voxel grid model to obtain the local neutron dose corresponding to each thread; Step S34, determine whether there are unselected source particles; Step S35: If there are unselected source particles, select a set number of source particles from the unselected source particles, and execute Steps S32 to S35 until all source particles are simulated in the target voxel grid model. Based on the multiple local neutron doses corresponding to each thread, determine the neutron dose distribution result of the target voxel grid model.

5. A method for calculating neutron dose according to claim 4, characterized in that, Based on multiple threads of the GPU itself, perform parallel simulation of all source particles in the batch of source particles to be simulated in the target voxel grid model, and obtain the local neutron dose corresponding to each thread, including: Based on a preset particle distribution rule, allocate each source particle in the batch of source particles to be simulated to the threads of the GPU itself, where each thread is allocated a selected source particle; Based on a preset simulation process, perform parallel simulation of all source particles in the batch of source particles to be simulated in the target voxel grid model, and calculate the energy deposition of each source particle in each voxel to obtain the local neutron dose corresponding to each thread.

6. A method for calculating neutron dose according to claim 5, characterized in that, The parallel simulation of all source particles in the batch of source particles to be simulated in the target voxel grid model based on the preset simulation process includes: Step S331: For each thread, calculate the reaction cross-section of the source particle corresponding to the thread in the target voxel model, and the value of the reaction cross-section characterizes the probability of the source particle undergoing a nuclear reaction with the target nucleus in the target voxel model; Step S332: For each thread, based on the transport length sampling kernel and the reaction cross-section of the source particle corresponding to the thread, perform transport length sampling of the source particle corresponding to the thread in the target voxel model to obtain the transport length of the source particle, and the value of the transport length characterizes the distance that the source particle moves in the target voxel model; Step S333: For each thread, based on the boundary distance calculation kernel, calculate the nearest boundary distance of the source particle corresponding to the thread along the movement direction to the voxel where it is located; Step S334: For each thread, based on the absorption reaction kernel, perform absorption processing on the source particle corresponding to the thread in the target voxel model to obtain the absorption state of the source particle corresponding to the thread; Step S335: For each thread, based on the scattering reaction kernel, perform scattering processing on the source particle corresponding to the thread in the target voxel model to obtain the scattering information of the source particle corresponding to the thread; Step S336: For each thread, based on the KERMA calculation kernel, calculate the KERMA values of multiple elements of the source particle corresponding to the thread in the voxel where it is located, and the KERMA value characterizes the initial kinetic energy of the charged particles ionized per unit mass of the target material by the source particle; Step S337: For each thread, based on the dose counting kernel and the KERMA value of the source particle corresponding to the thread, perform dose calculation on the source particle corresponding to the thread to obtain the local neutron dose corresponding to the thread. Step S338: For each of the threads, based on the transport length, the distance to the nearest boundary, the absorption state, and the scattering information of the source particles corresponding to the thread, determine whether there are source particles in the target voxel model; Step S339: For each of the threads, if there are still source particles in the target voxel model, execute Steps S332 to S339 until there are no particles in the target voxel model, and the simulation process ends.

7. Use of the method for calculating neutron dose according to any one of claims 1-6 in boron neutron capture therapy.

8. A neutron dose calculation device, characterized in that, Comprising: An acquisition module, configured to acquire an initial voxel grid model generated based on medical image information, the initial voxel grid model including a plurality of voxels, the plurality of voxels including air voxels, and the air voxels representing the volume elements occupied by air; A cutting module, configured to cut all the air voxels in the initial voxel grid model to generate a target voxel grid model; A simulation module, configured to acquire a plurality of source particles for Monte Carlo simulation, and based on the event-based Monte Carlo algorithm, simulate the plurality of source particles in the target voxel grid model to obtain the neutron dose distribution result of the target voxel grid model, and the neutron dose distribution result includes the neutron dose of each voxel in the target voxel grid model.

9. An electronic device, characterized in that, Comprising a processor, the processor being coupled to a memory; The processor is configured to execute a computer program stored in the memory to cause the electronic device to execute the method according to any one of claims 1 to 6.

10. A computer-readable storage medium, characterized in that, Comprising a computer program or instruction, when the computer program or instruction runs on a computer, causing the computer to execute the method according to any one of claims 1-6.

Images