Dose calculation method in boron neutron capture therapy planning system based on track tracking Monte Carlo method

Abstract

A method for dose calculation in a boron neutron capture therapy planning system based on a track tracking Monte Carlo method comprises the following steps: sampling according to input description source item particles, carrying out transport calculation by using the track tracking Monte Carlo method under voxel grid geometry, calculating a transport distance while calculating the transport distance, and calculating the dose in the boron neutron capture therapy planning system by using the track tracking Monte Carlo method. Counting based on the voxel grid to obtain particle flux and dose; by adopting the operation of track tracking method acceleration based on voxel grid geometry, the surface crossing judgment method is improved, and the problem that the surface crossing operation takes too long time in the calculation process is solved; counting is carried out while the transportation distance is calculated, the counting time is saved on the premise that the counting precision is guaranteed, and compared with a traditional Monte Carlo method based on CSG geometry, the method is more suitable for a voxel model obtained from a medical image, and the calculation speed is higher.

Description

Technical Field

[0001] The present invention relates to the field of treatment planning system (TPS) and particle transport calculation technology applied to boron neutron capture therapy (BNCT), and in particular to a method for dose calculation in a boron neutron capture therapy planning system based on a track tracing Monte Carlo method. Background Art

[0002] Boron neutron capture therapy (BNCT) is a new method for treating cancer. The method is to adjust the neutron flux generated by the accelerator or reactor neutron source to a specific energy range and spatial range through a beam shaper, irradiate the B-10 accumulated inside the tumor cells, induce the capture reaction of neutrons and B-10, and produce high-energy and high linear energy transfer (LET) α particles and Li-7 particles. The generated particles deposit a large amount of energy in a small range, thereby targeting and killing cancer cells under the premise of low deposition doses in normal cells. However, during the transport process before the neutrons are captured by B-10, doses will also be deposited in normal tissues. Therefore, it is necessary to predict the dose deposited by the treatment plan before treatment to optimize the treatment plan, so as to ensure that the biological equivalent dose in the tumor tissue is as large as possible under the premise that the dose to normal tissue does not exceed the limit. The system that predicts dose distribution and optimizes the treatment plan is called the Boron Neutron Capture Treatment Planning System (TPS), which mainly uses the Monte Carlo method or the deterministic method for dose prediction.

[0003] Due to the complexity of the reaction between neutrons and human tissues, the deterministic method can only simplify the reaction to make the equation solvable. Therefore, the calculation accuracy of the deterministic method is low and it is not widely used at present. The Monte Carlo method statistically samples particles, simulates their transport process in the human body, and obtains the particle flux in the counting grid. The obtained flux is multiplied by the flux-dose conversion factor to obtain the dose distribution of tissues in the human body. Since the Monte Carlo method is a statistical method, it can be determined according to the law of large numbers and the central limit law in mathematics that when the sample size is large enough, it can converge to the true solution.

[0004] When using the Monte Carlo method to predict BNCT dose distribution, particle transport is required based on the human body model. In the traditional Monte Carlo method, surface penetration judgment consumes a lot of time, accounting for about 70% of the total transport time. When constructing a geometric model for BNCT, the geometric and material information obtained by radiological medical imaging methods such as CT, MRI, and PET is based on pixels, so the human body model constructed based on this is a geometric model of structured voxels. In the traditional Monte Carlo method based on CSG geometry, this model requires a large number of surface penetration judgments for particles, resulting in low computational efficiency. Therefore, it is necessary to introduce a dedicated track tracking method to improve computational efficiency. Summary of the invention

[0005] In order to overcome the above problem, that is, the time-consuming particle penetration judgment in the prior art, resulting in low efficiency, the present invention aims to provide a method for dose calculation in a boron neutron capture therapy planning system based on a track tracking Monte Carlo method, which is based on the structural voxel grid geometry and uses the track tracking method to accelerate particle transport, thereby solving the problem of time-consuming particle penetration calculation.

[0006] In order to achieve the above object, the present invention adopts the following technical solutions to implement:

[0007] A method for dose calculation in a boron neutron capture therapy planning system based on a track-tracing Monte Carlo method, the steps are as follows:

[0008] Step 1: Establish a voxel-level counting grid based on the medical image, hereinafter referred to as the voxel grid, and set the counting result to zero;

[0009] Step 2: Sample source particles according to the source distribution parameters defined by the input card to obtain the position, energy and movement direction of the source particles;

[0010] Step 3: Determine the voxel grid position of the current particle according to the initial coordinates of the particle, and determine the material of the particle in the grid according to the grid position;

[0011] Step 4: According to the particle type, particle energy and the material in the grid where the particle is located, the particle flux dose conversion coefficient and macroscopic cross section are calculated, and the collision distance is sampled according to the macroscopic cross section;

[0012] Step 5: Based on the particle's position and direction of movement, calculate the distance from the particle's current position to the x-plane, y-plane, and z-plane of the next voxel grid;

[0013] Step 6: According to the direction of particle movement, calculate the particle's movement from the grid along the current transport direction. x face, y face or z The face crosses an entire voxel grid and reaches the next voxel grid x face, y face, or z The distance the surface is transported;

[0014] Step 7: Perform a through-surface iterative calculation to obtain the particle track length in the voxel grid, and determine whether the current voxel grid is a counting grid. If so, multiply the particle track length in the voxel grid by the particle flux dose conversion coefficient calculated in step 4 to obtain the deposition dose in each voxel grid, and accumulate the deposition dose in the voxel grid;

[0015] Step 8: Determine whether the cross-surface iteration is finished. If not, return to step 7. If the end condition is that the particle is absorbed by collision or transported out of the voxel grid boundary, go to step 9. If the end condition is that the material in the grid where the particle is located changes, or the particle changes its energy due to collision, return to step 4.

[0016] Step 9: Check whether there are any remaining particles in the secondary particle library. If yes, go back to step 3, otherwise go to step 10;

[0017] Step 10: Check whether the number of source particles in the voxel grid has reached the convergence limit. If not, execute step 2. Otherwise, normalize the counting results in the voxel grid one by one to give the dose distribution.

[0018] Compared with the prior art, the present invention has the following outstanding advantages:

[0019] 1. Since the traditional method for judging the direction of penetration needs to make multiple judgments in one penetration behavior, and the penetration judgment and distance statistics need to operate on different variables, the present invention improves the method of penetration judgment by comparing the distances between particles and different surfaces. Compared with the traditional Monte Carlo method, it saves the time used for geometric penetration judgment, reduces the calculation time required for dose prediction in the treatment plan, and improves the calculation efficiency;

[0020] 2. Counting is performed while making surface penetration judgment. Compared with the traditional Monte Carlo method, it reduces the calculation time used for counting and improves the calculation efficiency while ensuring the accuracy of dose calculation.

[0021] In summary, the method of the present invention reduces the time taken by the boron neutron capture therapy planning system to perform dose calculation while ensuring that the dose calculation accuracy remains unchanged, so that the time taken for dose calculation is more in line with clinical requirements. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 The figure is a flow chart of the method of the present invention.

[0023] Figure 2 This is a graphical representation of the quantities calculated in steps 5 and 6.

[0024] Figure 3 This is the iterative flow chart for surface penetration. DETAILED DESCRIPTION

[0025] The present invention is further described in detail below through the accompanying drawings and specific embodiments.

[0026] The present invention is a method for dose calculation in a boron neutron capture therapy planning system based on a track-tracing Monte Carlo method. Figure 1 As shown, the following technical solutions are adopted for implementation:

[0027] Step 1: Establish a voxel-level counting grid based on the medical image, that is, construct a set of counting results corresponding to the voxels in the area to be counted, and set the counting results to zero.

[0028] tally sum (n txmin :n txmax ,n tymin :n tymax ,n tzmin :n tzmax )=0 formula (1)

[0029] Where:

[0030] tally sum : is the total count amount in each cell, the length start coordinate and end coordinate of the three dimensions are the grid ordinal numbers corresponding to the start counting grid and the end counting grid in the voxel grid respectively;

[0031] n txmin , n tymin , n tzmin : are the grid numbers corresponding to the counting grids in the x, y, and z voxel grids at the beginning;

[0032] n txmax , n tymax , n tzmax : The grid numbers corresponding to the counted grids in the x, y, and z voxel grids at the end;

[0033] Step 2: Sample source particles according to the source distribution parameters defined by the input card to obtain the position, energy and movement direction of the source particles;

[0034] Step 3: Determine the voxel grid position of the current particle based on the initial coordinates of the particle (x0, y0, z0), that is, (n x ,n y ,n z ), and determine the material in the grid where the particle is located according to the grid position;

[0035] Step 4: According to the particle type, particle energy and the material in the grid where the particle is located, calculate the particle flux dose conversion coefficient K(P, E, M) and macroscopic cross section Σ(P, E, M), and sample the collision distance l according to the macroscopic cross section. c ;

[0036] Step 5: Based on the particle's position and direction of movement, calculate the distance from the particle's current position to the x-plane, y-plane, and z-plane of the next voxel grid in the voxel grid, as l x , l y and l z The initial value of lx0 , l y0 , l z0 , the initial value in two-dimensional geometry is Figure 2 As shown, A(x0,y0) is the initial position of the particle, and the particle moves to the upper right. x0 is the length of line segment AB1, l y0 is the length of line segment AB2; the situation is similar in three-dimensional geometry, the specific algorithm is as shown in formula (1), and the transported distance and the last transported distance l are recorded at the same time trans =0, l0=0:

[0037]

[0038] Where:

[0039] l x , l y and l z — Along the particle transport direction, greater than the transported distance l trans The initial value is the distance from the particle's starting position to the nearest voxel grid x-plane, y-plane, and z-plane. The initial value is the distance from the particle's initial position to the nearest voxel grid x-plane, y-plane, and z-plane.

[0040] (x0, y0, z0)—the initial position coordinates of the particle;

[0041] x t ,y t 、z t —The nearest voxel grid interface reached by the particle along the particle transport direction x Surface coordinates, y Surface coordinates and z Surface coordinates;

[0042] k x , k y , k z —The cosine value of the angle between the particle transport direction and the x-axis, y-axis, and z-axis;

[0043] Step 6: Find the particle's displacement from the grid along the current transport direction. x face, y face or z The face crosses an entire voxel grid and reaches the next voxel grid x face, y face, or z The distance d transported by the surface x ,d y or z , in two-dimensional geometry Figure 2 As shown, d x is the length of line segment B1C1, d yis the length of line segment B2C2; the situation is similar in three-dimensional geometry, and the specific algorithm is shown in formula (2):

[0044]

[0045] Where:

[0046] d x ,d y ,d z —From the voxel grid x noodle, y noodle, z The distance that the surface travels across an x-grid, y-grid, and z-grid along the particle transport direction;

[0047] Δx, Δy, Δz—the distance between two voxel grid surfaces;

[0048] k x , k y , k z —Particle transport direction and x axis, y axis, z The cosine of the axis angle;

[0049] Step 7: Record l x , l y and l z The minimum value in the index is index, and the transported distance at this time is l trans = l index , if the voxel grid (n x ,n y ,n z ) is the counting grid, and l tally =(l trans -l0)×K(P,E,M), that is, the dose generated by the particle in the grid during this transport, accumulated in the voxel grid (n x ,n y ,n z ) trans Assign it to l0, and then index To update, the update method is to add the corresponding d index , that is, according to the selected index, give l x , l y or z Add the corresponding d x ,d y or z ; At the same time, n index To update, the update method is to give n index According to the selected index, the voxel grid (n x ,ny ,n z ) x 、n y or z According to k x , k y or k z The positive or negative value is plus or minus 1.

[0050] Step 8: Return to step 7 to update l trans , until one of the following three situations occurs:

[0051] Case 1: Particle transport reaches the sampling collision distance, i.e. l trans ≥l c , at this time, the transport distance l trans Assigned to l c , and then count, that is, the particle transports l c After reaching a distance of , the particle collides with the material in the grid where the particle is located. If secondary particles are generated, the secondary particles are stored in the secondary particle library, and the movement direction and particle energy are changed according to the collision conditions, and the process returns to step 4. Or, due to collision truncation, the particle transport is terminated and the process goes to step 9.

[0052] Case 2: Particles are transported to the voxel grid boundary, that is, the updated voxel grid (n x ,n y ,n z ) is less than the minimum value when the grid is defined or greater than the maximum value when the grid is defined, the particle is transported out of the voxel grid, the particle transport ends, and the process goes to step 9;

[0053] Case 3: Particle transport to other materials, i.e. updated voxel grid (n x ,n y ,n z ) is different from the initial material, return to step 4, update the macroscopic cross section according to the material, and obtain a new collision distance;

[0054] The flowchart of the process is as follows Figure 3 shown.

[0055] Step 9: Check whether there are secondary particles in the secondary particle library. If so, take one out and then execute step 3 for the particle. If there are no secondary particles, the transport process of a single source particle is completed.

[0056] Step 10: According to the source term definition, return to step 2, re-extract particles, and execute steps 3 to 9 until the number of source term particles reaches the required given convergence limit, thereby ensuring the convergence of the calculation results. Normalize the results in the voxel grid one by one, give the dose distribution, obtain the distribution of absolute dose, and end the calculation.

Claims

1. A method for dose calculation in a boron neutron capture therapy planning system based on a track-tracing Monte Carlo method, characterized in that: Here are the steps: Step 1: Establish a voxel-level counting grid based on the medical image, hereinafter referred to as the voxel grid, and set the counting result to zero; Step 2: Sample source particles according to the source distribution parameters defined by the input card to obtain the position, energy and movement direction of the source particles; Step 3: Determine the voxel grid position of the current particle according to the initial coordinates of the particle, and determine the material of the particle in the grid according to the grid position; Step 4: According to the particle type, particle energy and the material in the grid where the particle is located, the particle flux dose conversion coefficient and macroscopic cross section are calculated, and the collision distance is sampled according to the macroscopic cross section; Step 5: Based on the particle's position and direction of movement, calculate the distance from the particle's current position to the x-plane, y-plane, and z-plane of the next voxel grid; Step 6: According to the direction of particle movement, calculate the particle's movement from the grid along the current transport direction. x face, y face or z The face crosses an entire voxel grid and reaches the next voxel grid x face, y face, or z The distance the surface is transported; Step 7: Perform a through-surface iterative calculation to obtain the particle track length in the voxel grid, and determine whether the current voxel grid is a counting grid. If so, multiply the particle track length in the voxel grid by the particle flux dose conversion coefficient calculated in step 4 to obtain the deposition dose in each voxel grid, and add the deposition dose to the voxel grid established in step 1; Step 8: Determine whether the cross-surface iteration is finished. If not, return to step 7. If the end condition is that the particle is absorbed by collision or transported out of the voxel grid boundary, go to step 9. If the end condition is that the material in the grid where the particle is located changes, or the particle changes its energy due to collision, return to step 4. Step 9: Check whether there are any remaining particles in the secondary particle library. If yes, go back to step 3, otherwise go to step 10; Step 10: Check whether the number of source particles in the voxel grid has reached the convergence limit. If not, execute step 2. Otherwise, normalize the counting results in the voxel grid one by one to give the dose distribution.

Images