A method for calculating proton radiation dose in a magnetic field

Abstract

A method for calculating proton radiation dose in a magnetic field comprises: determining the position of an emission source and the velocity and energy of a pencil beam at the emission source; after the proton beam is emitted from the emission source, determining multiple candidate boxes in a time sequence based on the position of the emission source, the velocity and energy of the pencil beam at the emission source, and information about the magnetic field, and sequentially calculating the position of the intersection of the pencil beam with the surface of each candidate box and the velocity and energy of the pencil beam at each intersection; selecting multiple intersection points in the candidate boxes inside a phantom as sampling points, connecting all the sampling points to obtain an actual motion trajectory of the pencil beam in the phantom; and obtaining the proton dose delivered by the pencil beam to any voxel in the phantom based on the actual motion trajectory and a pencil beam algorithm, thereby obtaining the proton dose irradiated to any voxel in the phantom.

Description

Technical Field

[0001] The present invention relates to the field of nuclear medicine imaging technology, and in particular to a method for calculating proton radiation dose in a magnetic field. Background Art

[0002] Proton therapy is a novel radiation therapy method for cancer treatment that uses high-energy proton beams. It offers advantages such as a relatively low input dose and sharp dose falloff, enabling precise dose delivery to the target while minimizing damage to surrounding normal tissue. However, factors such as proton range uncertainty, positioning errors, and tissue motion require physicists to impose margins around the target volume during radiotherapy planning, limiting the high dose gradients offered by proton therapy.

[0003] Magnetic resonance imaging (MRI) has advantages such as real-time imaging and high soft tissue contrast, and it does not emit ionizing radiation, making it particularly suitable for radiotherapy guidance of moving soft tissue tumors. Through real-time MRI guidance, the targeting accuracy of proton therapy, especially for moving soft tissue tumors, can be improved. However, to date, MRI-guided proton therapy (MRgPT) has not been achieved clinically. During the dose calculation, optimization, and delivery process, the dose distortion of the proton beam caused by the magnetic field effect (i.e., the Lorentz force) is one of the many obstacles to achieving MRgPT.

[0004] Wolf and Bortfeld first derived an analytical solution for proton beam deflection trajectories in a transverse uniform magnetic field. Their method is based on a small-angle approximation, resulting in increased errors for large proton deflection angles; furthermore, their method lacks a closed-form solution in the relativistic case. Schellhammer and Hoffmann introduced an analytical iterative method that yields trajectories in good agreement with Monte Carlo simulations in the energy range of 60 to 250 MeV. However, the method's applicability is limited to cases where the magnetic field is perpendicular to the proton beam and is not applicable to voxelized phantoms. Padilla-Cabal proposed a numerical iterative method that stores the energy deposition values ​​of the proton beam in water along the trajectory in a lookup table. For non-aqueous materials, this method determines the water-equivalent depth using a lookup table and scales it using a material-specific correction factor. However, this correction factor does not accurately describe the properties of high-density materials such as bone, leading to biased dose distributions in these materials. Monte Carlo simulations provide high-precision simulations, but their computational speed limits their application in real-time MRgPT. Therefore, an efficient and relatively accurate dose calculation algorithm is needed to better realize the radiation dose calculation in MRgPT. Summary of the Invention

[0005] In view of the above problems, the present invention provides a method for calculating proton radiation dose in a magnetic field to achieve efficient and relatively accurate determination of proton radiation dose.

[0006] As a first aspect of the present invention, a method for calculating a proton radiation dose in a magnetic field is provided, wherein a proton beam is emitted from an emission source and enters the magnetic field, the proton beam moves in the magnetic field and enters a phantom located in the magnetic field to irradiate protons into the phantom, the proton beam comprising a plurality of pencil beams, the method comprising:

[0007] determining a position of the emission source and a velocity and energy of the pencil beam at the emission source;

[0008] After the proton beam is emitted from the emission source, a plurality of candidate boxes are determined in a time sequence according to the position of the emission source, the velocity and energy of the pencil beam at the emission source, and information about the magnetic field, and the position of the intersection of the pencil beam and the surface of each candidate box and the velocity and energy of the pencil beam at each intersection are calculated in sequence, and the candidate box corresponding to the point where the intersection is outside the phantom or the energy at the intersection is zero is selected as the last candidate box. Of two adjacent candidate boxes, the later candidate box is determined according to the position of the intersection of the pencil beam and the earlier candidate box and the velocity of the pencil beam at the intersection.

[0009] In the candidate box inside the phantom, multiple intersection points are selected as sampling points, and all the sampling points are connected to obtain the actual motion trajectory of the pencil beam in the phantom;

[0010] The proton dose delivered to any voxel in the phantom by the pencil beam is obtained according to the actual motion trajectory and the pencil beam algorithm, and the proton dose irradiated to any voxel in the phantom by the proton beam is further obtained.

[0011] According to an embodiment of the present invention, after the proton beam is emitted from the emission source, a plurality of candidate boxes are determined in a time sequence according to the position of the emission source, the velocity and energy of the pencil beam at the emission source, and information of the magnetic field, and the position of the intersection point of the pencil beam with the surface of each candidate box and the velocity and energy of the pencil beam at each intersection point are calculated in sequence, including:

[0012] Use {i|0≤i≤H} to number multiple candidate boxes and multiple intersection points determined in time sequence; use the phantom as the 0th candidate box and the position of the emission source as the 0th intersection point;

[0013] An i-th candidate box is determined based on the CT image data of the phantom, the position of the i-th intersection, the velocity and energy of the pencil beam at the i-th intersection, and the information of the magnetic field, and the position of the i+1-th intersection of the pencil beam and the i-th candidate box and the energy and velocity of the pencil beam at the i+1-th intersection are calculated. When i=H, the H-th intersection is outside the phantom, or the velocity of the pencil beam at the H-th intersection is zero, and the i+1-th intersection is not calculated.

[0014] According to an embodiment of the present invention, the angle between the pencil beam velocity direction and the magnetic field direction is θ, and θ∈[0,π];

[0015] A method for calculating a position of an (i+1)th intersection point where the pencil beam intersects the surface of the (i)th candidate box, comprising:

[0016] determining a theoretical motion trajectory of the pencil beam from the i-th intersection point according to the position of the i-th intersection point, the velocity of the pencil beam at the i-th intersection point, and information of the magnetic field;

[0017] The position of the (i+1)th intersection point is obtained according to a theoretical motion trajectory of the pencil beam from the (i)th intersection point and a boundary condition of the (i)th candidate box.

[0018] According to an embodiment of the present invention, the theoretical motion trajectory is a spiral trajectory.

[0019] Determining a theoretical motion trajectory of the pencil beam from the i-th intersection point according to the position of the i-th intersection point, the velocity of the pencil beam at the i-th intersection point, and information of the magnetic field, comprising:

[0020] Obtaining an initial phase angle of the spiral trajectory of the pencil beam starting from the i-th intersection point according to the position of the i-th intersection point;

[0021] determining a radius and a pitch of a spiral trajectory of the pencil beam starting from the i-th intersection point according to the velocity of the pencil beam at the i-th intersection point, information about the magnetic field direction, and the energy of the pencil beam at the i-th intersection point;

[0022] The theoretical motion trajectory of the pencil beam starting from the i-th intersection is obtained according to the initial phase angle, radius and pitch of the spiral trajectory of the pencil beam starting from the i-th intersection.

[0023] According to an embodiment of the present invention, at each intersection point in the magnetic field, the velocity of the pencil beam is corrected using a coordinate basis matrix.

[0024] According to an embodiment of the present invention, a method for determining the energy and velocity of the pencil beam at the (i+1)th intersection point includes:

[0025] Obtaining Bragg-Kleeman parameters of the i-th candidate box in the phantom according to the CT image data of the phantom;

[0026] The energy and speed of the pencil beam at the (i+1)th intersection are determined according to the Bragg-Kleeman parameters of the (i)th candidate box, the speed and energy of the pencil beam at the (i)th intersection.

[0027] According to an embodiment of the present invention, obtaining the Bragg-Kleeman parameters of the i-th candidate box in the phantom according to the CT image data of the phantom includes:

[0028] Make the voxel of the CT image the same size as the i-th candidate box, where i>0;

[0029] Extracting the Bragg-Kleeman parameter of the i-th candidate box from a Bragg-Kleeman parameter lookup table according to the HU value of the voxel in the CT image;

[0030] According to an embodiment of the present invention, the magnetic field is a uniform magnetic field, and obtaining a theoretical motion trajectory of the pencil beam from the i-th intersection point based on the initial phase angle, radius, and pitch of the spiral trajectory of the pencil beam from the i-th intersection point includes:

[0031] Establishing a spiral parameter equation of the pencil beam starting from the i-th intersection in a local coordinate system using the initial phase angle, radius, and pitch of the spiral trajectory of the pencil beam starting from the i-th intersection;

[0032] The spiral parameter equation of the pencil beam from the i-th intersection is converted to a global coordinate system to obtain the spiral parameter equation of the pencil beam from the i-th intersection in the global coordinate system. The spiral parameter equation of the pencil beam from the i-th intersection in the global coordinate system is used to describe the theoretical motion trajectory of the pencil beam from the i-th intersection.

[0033] According to an embodiment of the present invention, obtaining the position of the (i+1)th intersection point according to the theoretical motion trajectory of the pencil beam from the (i)th intersection point and the boundary condition of the (i)th candidate box includes:

[0034] Taking three planes in the i-th candidate box that have a probability of intersecting with the pencil beam as candidate planes;

[0035] Substituting the boundary conditions of the three candidate planes into the spiral parametric equation of the pencil beam starting from the i-th intersection in the global coordinate system, a unique valid solution is obtained;

[0036] Substitute the effective solution back into the spiral parameter equation of the pencil beam starting from the i-th intersection in the global coordinate system to obtain the position of the i+1-th intersection.

[0037] According to an embodiment of the present invention, the geometric center of the phantom coincides with the geometric center of the magnetic field.

[0038] According to embodiments of the present invention, by setting up multiple candidate boxes and calculating the intersection of the pencil beam with the surface of each candidate box, the actual trajectory of the pencil beam in a phantom is determined. Combining the actual trajectory with the pencil beam algorithm, the proton dose delivered by the pencil beam to any voxel in the phantom is calculated, thereby enabling determination of the proton radiation dose in a magnetic field. The method provided by embodiments of the present invention has high accuracy and can be used to implement magnetic resonance-guided proton therapy for patients. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 A method for calculating proton radiation dose in a magnetic field according to an embodiment of the present invention is shown;

[0040] Figure 2 A schematic diagram of an intersection provided according to an embodiment of the present invention is shown;

[0041] Figure 3 A schematic diagram showing a theoretical motion trajectory of a pencil beam in a magnetic field in a local coordinate system according to an embodiment of the present invention is shown;

[0042] Figure 4 The dose distribution of the proton beam incident on the prostate with a gantry angle of 68° / 248° under a 3.0T uniform transverse magnetic field.

[0043] Figure 5 Dose distribution of a proton beam incident on a water phantom at a 3.0 T uniform transverse magnetic field and a gantry angle of 270°. DETAILED DESCRIPTION

[0044] In order to make the objectives, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with specific embodiments and with reference to the accompanying drawings. Obviously, the embodiments described are part of the embodiments of the present invention, not all of them. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0045] Figure 1 A method for calculating proton radiation dose in a magnetic field according to an embodiment of the present invention is shown. A proton beam is emitted from a source and enters the magnetic field. The proton beam moves within the magnetic field and enters a phantom within the magnetic field to irradiate protons toward the phantom. The proton beam includes multiple pencil beams. The calculation method includes steps S1 to S4.

[0046] In step S1 , the position of the emission source and the speed and energy of the pencil beam at the emission source are determined.

[0047] In step S2, after the proton beam is emitted from the emission source, multiple candidate boxes are determined in a time sequence based on the position of the emission source, the velocity and energy of the pencil beam at the emission source, and information about the magnetic field. The position of the intersection of the pencil beam and the surface of each candidate box and the velocity and energy of the pencil beam at each intersection are calculated in sequence. The candidate box corresponding to the point where the intersection is outside the phantom or the energy at the intersection is zero is selected as the last candidate box. Among two adjacent candidate boxes, the later candidate box is determined based on the position of the intersection of the pencil beam and the earlier candidate box and the velocity of the pencil beam at the intersection.

[0048] In step S3, multiple intersection points are selected as sampling points in the candidate box inside the phantom, and all sampling points are connected to obtain the actual motion trajectory of the pencil beam in the phantom;

[0049] In step S4, the proton dose delivered by the pencil beam to any voxel in the phantom is obtained according to the actual motion trajectory and the pencil beam algorithm, and then the proton dose irradiated by the proton beam to any voxel in the phantom is obtained.

[0050] According to embodiments of the present invention, by setting up multiple candidate boxes and calculating the intersection of the pencil beam with the surface of each candidate box, the actual trajectory of the pencil beam in a phantom is determined. Combining the actual trajectory with the pencil beam algorithm, the proton dose delivered by the pencil beam to any voxel in the phantom is calculated, thereby enabling determination of the proton radiation dose in a magnetic field. The method provided by embodiments of the present invention has high accuracy and can be used to implement magnetic resonance-guided proton therapy for patients.

[0051] According to an embodiment of the present invention, in step S2, after the proton beam is emitted from the emission source, multiple candidate boxes are determined in a time sequence based on the position of the emission source, the velocity and energy of the pencil beam at the emission source, and information about the magnetic field, and the position of the intersection of the pencil beam with the surface of each candidate box and the velocity and energy of the pencil beam at each intersection are calculated in sequence, including steps S21 to S22.

[0052] In step S21, multiple candidate boxes and multiple intersection points determined in time sequence are numbered using {i|0≤i≤H}; the phantom is used as the 0th candidate box, and the position of the emission source is used as the 0th intersection point;

[0053] In step S22, an i-th candidate box is determined based on the CT image data of the phantom, the position of the i-th intersection point, the velocity and energy of the pencil beam at the i-th intersection point, and the magnetic field information. The position of the i+1-th intersection point of the pencil beam with the i-th candidate box and the energy and velocity of the pencil beam at the i+1-th intersection point are calculated. When i=H, the H-th intersection point is outside the phantom, or the velocity of the pencil beam at the H-th intersection point is zero, and the i+1-th intersection point is no longer calculated.

[0054] Figure 2 A schematic diagram of an intersection provided according to an embodiment of the present invention is shown.

[0055] like Figure 2 As shown, P i+1 represents the intersection of the pencil beam and the i-th candidate box, that is, the i+1-th intersection, and the direction of the arrow represents the velocity direction of the pencil beam at each intersection.

[0056] According to an embodiment of the present invention, in step S1, the emission source can be located either within or outside the magnetic field, depending on the magnitude of the magnetic field. The emission source in step S1 is not a physical point source within the treatment machine head. The proton beam is generated at a physical point source, passes through various components within the proton therapy machine head, and is scanned and diffused into a wider two-dimensional plane, which is the emission source. The emission source is located at the trailing edge of the range shifter and can also be referred to as a virtual source.

[0057] According to an embodiment of the present invention, the velocity of the pencil beam at the exit source includes magnitude and direction, and the velocity and energy of the pencil beam at the exit source and the position of the exit source can be obtained based on the radiotherapy planning file. The information provided by the radiotherapy planning file includes the number of irradiation points (i.e., the number of pencil beams covering the phantom), the gantry angle, head position, isocenter, source wheelbase, and shifter information of the treatment machine in the beam viewing coordinate system, the angle of the treatment bed used to place the phantom, and the irradiation time (in MU) of each control point (pencil beam) and the nominal energy of the pencil beam. The initial phase space information can be obtained based on the radiotherapy planning file. The phase space information includes the coordinates of the pencil beam in the International Electrotechnical Commission (IEC) coordinate system (i.e., the initial position of the pencil beam), the direction vector (i.e., the velocity direction of the pencil beam at the exit source), the divergence angle, the beam spot size, and the correlation coefficient. The energy of the pencil beam at the exit source can be obtained using the beam energy model.

[0058] According to an embodiment of the present invention, the angle between the pencil beam velocity direction and the magnetic field direction is θ, and θ∈[0,π]. In step S22, a method for calculating the position of the (i+1)th intersection point where the pencil beam intersects the surface of the (i)th candidate box includes steps S221 and S222.

[0059] Step S221, determining a theoretical motion trajectory of the pencil beam starting from the i-th intersection point based on the position of the i-th intersection point, the velocity of the pencil beam at the i-th intersection point, and the magnetic field information;

[0060] Step S222 : obtaining the position of the (i+1)th intersection point according to the theoretical motion trajectory of the pencil beam from the (i)th intersection point and the boundary condition of the (i)th candidate box.

[0061] According to an embodiment of the present invention, the theoretical motion trajectory is a spiral trajectory. In step S222,

[0062] According to the position of the i-th intersection, the velocity of the pencil beam at the i-th intersection, and the magnetic field information, a theoretical motion trajectory of the pencil beam starting from the i-th intersection is determined, including S2211-S2213.

[0063] Step S2211, obtaining the initial phase angle of the spiral trajectory of the pencil beam starting from the i-th intersection point according to the position of the i-th intersection point;

[0064] Step S2212, determining the radius and pitch of the spiral trajectory of the pencil beam starting from the i-th intersection point according to the velocity of the pencil beam at the i-th intersection point, the information of the magnetic field direction, and the energy of the pencil beam at the i-th intersection point;

[0065] Step S2213 , obtaining a theoretical motion trajectory of the pencil beam starting from the i-th intersection point according to the initial phase angle, radius, and pitch of the spiral trajectory of the pencil beam starting from the i-th intersection point.

[0066] According to an embodiment of the present invention, the magnetic field may be a uniform magnetic field or a non-uniform magnetic field. For example, in the case of a uniform magnetic field, the magnetic field strength may be 3.0 T, the magnetic field direction may be along the negative Z-axis of the IEC coordinate system, and the shape may be a rectangular parallelepiped (the length, width, and height may all be 50 cm). S2213 obtains a theoretical motion trajectory of the pencil beam from the i-th intersection point based on the initial phase angle, radius, and pitch of the spiral trajectory of the pencil beam from the i-th intersection point, including the following steps S22131-S22132.

[0067] Step S22131 : establishing a spiral parameter equation of the pencil beam starting from the i-th intersection in the local coordinate system using the initial phase angle, radius, and pitch of the spiral trajectory of the pencil beam starting from the i-th intersection.

[0068] According to an embodiment of the present invention, without considering the energy loss for the time being, the spiral parameter equation P is established in the local coordinate system (LCS) lcs (φ)=(x lcs (φ),y lcs (φ),z lcs (φ)) is shown below.

[0069]

[0070] According to an embodiment of the present invention, the radius of the spiral trajectory of the pencil beam starting from the i-th intersection point is a i , the pitch is b i , the initial phase angle of the spiral trajectory is φ 0,i

[0071] Step S22132: Convert the spiral parametric equation of the pencil beam from the i-th intersection point to the global coordinate system to obtain the spiral parametric equation of the pencil beam from the i-th intersection point in the global coordinate system. The spiral parametric equation of the pencil beam from the i-th intersection point in the global coordinate system is used to describe the theoretical motion trajectory of the pencil beam from the i-th intersection point.

[0072] Using the orthogonal matrix M composed of LGS basis vectors i , P lcs (φ) is converted to GCS, and the trajectory equation of the proton in GCS is obtained P(φ). According to an embodiment of the present invention, the parametric equation of the spiral trajectory (Formula (1)) can only be expressed in a local coordinate system. The Z axis of the local coordinate system is antiparallel to the direction of the magnetic field, and the orthogonal basis matrix of the local coordinate system is denoted as M i

[0073] According to an embodiment of the present invention, in step S22 , the method for determining the energy and velocity of the pencil beam at the i-th intersection point includes operations S223 to S224 .

[0074] Step S223 : obtaining the Bragg-Kleeman parameters of the i-th candidate box in the phantom according to the CT image data of the phantom.

[0075] Step S224 : determining the energy and velocity of the pencil beam at the (i+1)th intersection point according to the Bragg-Kleeman parameters of the i-th candidate box and the velocity and energy of the pencil beam at the i-th intersection point.

[0076] According to an embodiment of the present invention, obtaining the Bragg-Kleeman parameters of the i-th candidate box in the phantom according to the CT image data of the phantom in step S223 includes steps S2231 to S2235.

[0077] Step S2231 : Make the voxel of the CT image the same size as the i-th candidate box, where i>0.

[0078] Step S2232 : extracting the Bragg-Kleeman parameters of the i-th candidate box from the Bragg-Kleeman parameter lookup table according to the HU values ​​of the voxels in the CT image.

[0079] According to an embodiment of the present invention, the curvature of the spiral motion trajectory of a pencil beam in a magnetic field is related to its energy. To maximize the accuracy of deflection prediction, an example embodiment of the present invention uses a comparison table of HU values, Bragg-Kleeman parameters (α, p), and material mass density to account for the impact of energy changes on the relative stopping power (RSP) of protons during the calculation process. After HU values ​​in the range of -1000-2995 are mapped to 25 materials according to a parameter conversion file, a Monte Carlo program is used to simulate the integrated depth dose curves (IDD) of the pencil beam in the energy range of 10-250 MeV and the energy interval of 10 MeV in the 25 materials to obtain the range R at each nominal energy E. E and R are then fitted to obtain the Bragg-Kleeman parameters for each material. Embodiments of the present invention can also choose to use HU-RSP calibration curves to calculate the beam trajectory.

[0080] In step S224, the energy and velocity of the pencil beam at the i+1 intersection are determined based on the Bragg-Kleeman parameters of the i-th candidate box, the velocity and energy of the pencil beam at the i-th intersection, which can be expressed as follows. According to an embodiment of the present invention, the medium will lose the energy of the proton beam, so the Bragg-Kleeman rule is used to correct the residual energy E of the pencil beam at the i+1 intersection. i+1 :

[0081]

[0082] According to an embodiment of the present invention, wherein α i ,p i is the Bragg-Kleeman parameter of the i-th candidate box, ρ i is the mass density of the i-th candidate box, l arc is the distance traveled by the pencil beam in the i-th candidate box.

[0083] According to an embodiment of the present invention, due to the influence of the Lorentz force, the proton beam will be deflected in the velocity direction perpendicular to the magnetic field direction, so it is necessary to use a rotation matrix R rotating around the Z axis. z Correct the velocity direction. If the proton is in the air outside the medium, only the direction is corrected. If it is in the voxel, the energy is corrected as well. The velocity direction correction formula is as follows:

[0084]

[0085] in, represents the unit vector of the pencil beam velocity at the i+1th intersection point, represents the component of the unit vector of the pencil beam velocity at the i-th intersection point parallel to the magnetic field, represents the component of the unit velocity vector of the pencil beam at the i-th intersection point in the direction perpendicular to the magnetic field, M i The transposed matrix, R z (φ i+1 ) T R z (φ i+1 ), φ i+1 is the angle at which the pencil beam is deflected when it reaches the i+1th intersection point.

[0086] According to an embodiment of the present invention, the geometric center of the phantom coincides with the geometric center of the magnetic field.

[0087] When the proton source is outside the magnetic field, the theoretical trajectory of the pencil beam before reaching the phantom is determined by the location of the proton source, the geometric boundaries of the magnetic field, the direction and strength of the magnetic field, and the energy and direction of the pencil beam. This is because the pencil beam moves in a straight line before entering the magnetic field and only begins to move in a curved line after entering the field due to the Lorentz force. Therefore, when the proton source is outside the magnetic field, the trajectory of the pencil beam before entering the magnetic field cannot be determined using a spiral trajectory.

[0088] According to an embodiment of the present invention, in step S222 , the position of the (i+1)th intersection point is obtained according to the theoretical motion trajectory of the pencil beam from the (i)th intersection point and the boundary condition of the (i)th candidate box, including steps S2221 to S2223 .

[0089] Step S2221: three planes in the i-th candidate box that have a probability of intersecting with the pencil beam are selected as candidate planes;

[0090] Step S2222, substituting the boundary conditions of the three candidate planes into the spiral parametric equation of the pencil beam starting from the i-th intersection in the global coordinate system to obtain a unique valid solution;

[0091] In step S2223 , the effective solution is substituted back into the spiral parameter equation of the pencil beam starting from the i-th intersection in the global coordinate system to obtain the position of the i+1-th intersection.

[0092] According to an embodiment of the present invention, the geometric body representing the phantom is taken as the 0th candidate box, and the three surfaces facing the proton momentum direction in the pencil beam are taken as candidate planes; the voxel i where the pencil beam is located is taken as the i-th candidate box, and the three surfaces facing away are taken as candidate planes.

[0093] According to an embodiment of the present invention, for the 0th candidate box, the solution interval is set to [0, 0.3π], and for the i-th candidate box (i>0), the solution interval is [0, the diagonal length of the i-th candidate box / a i ]; By judging the intersection point P i+1 Is it located on the surface of the candidate box? Get the only valid solution φ from the three solutionsi+1 .

[0094] According to an embodiment of the present invention, in step S3, the line connecting two adjacent sampling points among the n sampling points is a straight line, and the n sampling points are connected to form a broken line, and the phantom is approximately divided into n-1 calculation volumes. The number and spacing of the sampling points can be adjusted according to the requirements of calculation accuracy and efficiency.

[0095] According to an embodiment of the present invention, after step S3, the method further comprises calculating the physical path length l and radial distance r for each voxel in the phantom. The radial distance is the distance from the center of a voxel to the foot of a perpendicular line drawn through the actual motion trajectory. The physical path length is the distance from the first sampling point along the actual motion trajectory of the pencil beam to the foot of the perpendicular line. For an actual motion trajectory, each voxel has both a radial distance and a physical path length.

[0096] According to an embodiment of the present invention, in step S4, the proton motion trajectory is combined with the pencil beam algorithm, and the dose d(r,l,l) delivered by a pencil beam to a voxel located at the point of interest (r,l) is w ) can be calculated by multiplying the integrated depth dose (IDD) by the double Gaussian kernel function (K) using the following formula:

[0097] d(r,l,l w )= IDD(l w )K(r,l,l w ) (4)

[0098] K(r,l,l w )= (1-W nuc (l w ))G mcs (r,l)+W nuc (l w )G nuc (r,l w ) (5)

[0099]

[0100] Among them, l w is the equivalent water depth, W nuc is the nuclear reaction weight, G mcs is the contribution of multiple Coulomb scattering, G nuc is the contribution of nuclear reactions, σ mcs is the variance of the multiple Coulomb scattering Gaussian kernel, σ nuc is the variance of the Gaussian kernel of the nuclear reaction. mcs and σ nuc It is obtained through the divergence angle, beam spot size and correlation coefficient in the initial phase space information.

[0101] Among them, IDD(l w ) were extracted from the IDD database based on equivalent water depth. The IDD database covers a depth range of 0-40 cm with a resolution of 0.1 mm and an energy range of 10.1-250.0 MeV with an energy resolution of 0.1 MeV. This IDD database was constructed using Monte Carlo simulations.

[0102] According to an embodiment of the present invention, the doses of all pencil beams (PB) are superimposed to obtain a dose distribution D(x, y, z) of a voxel located at (x, y, z) under a magnetic field.

[0103]

[0104] Among them, w PB is the weight of each pencil beam.

[0105] According to an embodiment of the present invention, a Gaussian energy spectrum is derived from IDD measurement data of a pencil beam emitted by an actual treatment machine at different nominal energies. The beam's energy spread and central energy can be extracted from this Gaussian spectrum. The initial spot size, angular divergence, and correlation of the pencil beam are determined based on the treatment machine's air fluence (IAF) measurement data. The pencil beam is detected using a detector, and the number of emitted particles per unit (MU) of the pencil beam at different nominal energies is determined based on the dose deposition in the detector. This number of emitted particles per unit (MU) of the pencil beam is used to determine the weight of the pencil beam.

[0106] According to embodiments of the present invention, proton radiation dose determination in a magnetic field can also be performed based on MRI images. However, because MRI signal intensity depends on proton density and tissue relaxation properties, it cannot be directly used for dose calculation. MRI images must be converted into CT HU maps, also known as "pseudo-CT (sCT)" images. A mapping function between MRI voxel intensity and HU values ​​is constructed using extensive clinical data to achieve this conversion, allowing ion beam radiotherapy planning to be developed based on sCT.

[0107] Figure 4 Dose distribution of the proton beam incident on the prostate with a gantry angle of 68° / 248° under a 3.0T uniform transverse magnetic field.

[0108] Figure 4 Part (a) is the result determined by the embodiment of the present invention, part (b) is the simulation result of the Monte Carlo method, and part (c) is the 2mm / 2% Gamma pass rate diagram of the present invention. Figure 4 It can be seen that the consistency with the Monte Carlo simulation results demonstrates the accuracy of the present invention, and the 2mm / 2% Gamma passing rate at the 10% threshold is 99.20%.

[0109] Figure 5 Dose distribution of a proton beam incident on a water phantom at a 3.0 T uniform transverse magnetic field and a gantry angle of 270°.

[0110] Figure 5 Part (a) is the result determined by the embodiment of the present invention, part (b) is the simulation result of the Monte Carlo method, and part (c) is the 2mm / 2% Gamma pass rate diagram of the present invention. Figure 5 It can be seen that the consistency with the Monte Carlo simulation results demonstrates the accuracy of the present invention, and the 2mm / 2% Gamma passing rate at the 10% threshold is 98.97%.

[0111] In practical applications, such as when performing proton therapy on patients, it is necessary to obtain organ delineation files, patient CT images, and radiotherapy plan files.

[0112] The present invention utilizes local and global coordinate systems, based on the helical parametric equation and the relativistic Lorentz equation, to derive the theoretical trajectory equation for a pencil beam in a magnetic field. Substituting the candidate box surface boundary conditions, the trajectory of the pencil beam in a phantom under a magnetic field is numerically solved. Finally, combined with a pencil beam algorithm, the dose distribution of proton radiation under a magnetic field is calculated, thereby resolving three major challenges in this field. First, the applicability of existing ray tracing analytical algorithms for charged particles under a magnetic field is limited to cases where the magnetic field is perpendicular to the proton beam and is not applicable to voxelized phantoms. Second, existing numerical iterative algorithms for charged particles under a magnetic field cannot accurately describe the properties of high-density materials, resulting in large errors in the dose distribution in these materials. Third, existing traditional Monte Carlo programs have slow computational speeds and are unable to meet the speed requirements for real-time dose calculation in practical MRgRT. The present invention solves these problems and can be used to calculate proton radiation dose under a magnetic field in treatment planning systems, shortening calculation time while ensuring accuracy.

[0113] The method provided by the embodiment of the present invention is applicable to incident proton beams in any direction and any heterogeneous phantom, and can be applied to both uniform magnetic fields and non-uniform magnetic fields.

[0114] The method provided by the embodiment of the present invention only generates a single-variable nonlinear function when calculating the actual motion trajectory of the pencil beam. The solution range can be set according to the size of the candidate box, which greatly reduces the number of iterations required for the numerical solution and improves the calculation speed.

[0115] The method provided in the embodiment of the present invention samples the intersection points and approximates the proton trajectory as a broken line connecting the sampling points to calculate the physical path length l and radial distance r of each voxel in the medium, which is successfully combined with the pencil beam algorithm.

[0116] The above specific embodiments further illustrate the objectives, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above are only specific embodiments of the present invention and are not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. A method for calculating proton radiation dose in a magnetic field, wherein: A proton beam is emitted from an emission source and enters a magnetic field. The proton beam moves in the magnetic field and enters a phantom located in the magnetic field to irradiate protons to the phantom. The proton beam includes a plurality of pencil beams. The method includes: determining a position of the emission source and a velocity and energy of the pencil beam at the emission source; After the proton beam is emitted from the emission source, a plurality of candidate boxes are determined in a time sequence according to the position of the emission source, the velocity and energy of the pencil beam at the emission source, and information about the magnetic field, and the position of the intersection of the pencil beam and the surface of each candidate box and the velocity and energy of the pencil beam at each intersection are calculated in sequence, and the candidate box corresponding to the point where the intersection is outside the phantom or the energy of the intersection is zero is taken as the last candidate box, wherein, between two adjacent candidate boxes, the later candidate box is determined according to the position of the intersection of the pencil beam and the earlier candidate box and the velocity of the pencil beam at the intersection; In the candidate box inside the phantom, multiple intersection points are selected as sampling points, and all the sampling points are connected to obtain the actual motion trajectory of the pencil beam in the phantom; The proton dose delivered to any voxel in the phantom by the pencil beam is obtained according to the actual motion trajectory and the pencil beam algorithm, and the proton dose irradiated to any voxel in the phantom by the proton beam is further obtained.

2. The calculation method according to claim 1, wherein: After the proton beam is emitted from the emission source, a plurality of candidate boxes are determined in a time sequence according to the position of the emission source, the speed and energy of the pencil beam at the emission source, and information of the magnetic field, and the position of the intersection point of the pencil beam with the surface of each candidate box and the speed and energy of the pencil beam at each intersection point are calculated in sequence, including: Use {i|0≤i≤H} to number multiple candidate boxes and multiple intersection points determined in time sequence; use the phantom as the 0th candidate box and the position of the emission source as the 0th intersection point; An i-th candidate box is determined based on the CT image data of the phantom, the position of an i-th intersection, the velocity and energy of the pencil beam at the i-th intersection, and the information of the magnetic field, and the position of an i+1-th intersection of the pencil beam and the i-th candidate box and the energy and velocity of the pencil beam at the i+1-th intersection are calculated. When i=H, the H-th intersection is outside the phantom, or the velocity of the pencil beam at the H-th intersection is zero, and the i+1-th intersection is not calculated.

3. The calculation method according to claim 2, wherein: The angle between the pencil beam velocity direction and the magnetic field direction is ,and A method for calculating a position of an (i+1)th intersection point where the pencil beam intersects the surface of the (i)th candidate box, comprising: determining a theoretical motion trajectory of the pencil beam from the i-th intersection point according to the position of the i-th intersection point, the velocity of the pencil beam at the i-th intersection point, and information of the magnetic field; The position of the (i+1)th intersection point is obtained according to a theoretical motion trajectory of the pencil beam from the (i)th intersection point and a boundary condition of the (i)th candidate box.

4. The calculation method according to claim 3, wherein: The theoretical motion trajectory is a spiral trajectory. Determining a theoretical motion trajectory of the pencil beam from the i-th intersection point according to the position of the i-th intersection point, the velocity of the pencil beam at the i-th intersection point, and information of the magnetic field, comprising: Obtaining an initial phase angle of the spiral trajectory of the pencil beam starting from the i-th intersection point according to the position of the i-th intersection point; determining a radius and a pitch of a spiral trajectory of the pencil beam starting from the i-th intersection point according to the velocity of the pencil beam at the i-th intersection point, information about the magnetic field direction, and the energy of the pencil beam at the i-th intersection point; The theoretical motion trajectory of the pencil beam starting from the i-th intersection is obtained according to the initial phase angle, radius and pitch of the spiral trajectory of the pencil beam starting from the i-th intersection. 5 . The calculation method according to claim 2 , further comprising, at each intersection point in the magnetic field, correcting the velocity of the pencil beam using a coordinate basis matrix.

6. The calculation method according to claim 2, wherein: The method for determining the energy and velocity of the pencil beam at the (i+1)th intersection point comprises: Obtaining Bragg-Kleeman parameters of the i-th candidate box in the phantom according to the CT image data of the phantom; The energy and speed of the pencil beam at the (i+1)th intersection are determined according to the Bragg-Kleeman parameters of the (i)th candidate box, the speed and energy of the pencil beam at the (i)th intersection.

7. The calculation method according to claim 6, wherein: Obtaining the Bragg-Kleeman parameters of the i-th candidate box in the phantom according to the CT image data of the phantom includes: Make the voxel of the CT image the same size as the i-th candidate box, where i>0; The Bragg-Kleeman parameters of the i-th candidate box are extracted from a Bragg-Kleeman parameter lookup table according to the HU values ​​of the voxels in the CT image.

8. The calculation method according to claim 4, wherein: The magnetic field is a uniform magnetic field, and a theoretical motion trajectory of the pencil beam from the i-th intersection is obtained according to the initial phase angle, radius, and pitch of the spiral trajectory of the pencil beam from the i-th intersection, including: Establishing a spiral parameter equation of the pencil beam starting from the i-th intersection in a local coordinate system using the initial phase angle, radius, and pitch of the spiral trajectory of the pencil beam starting from the i-th intersection; The spiral parameter equation of the pencil beam from the i-th intersection is converted to a global coordinate system to obtain the spiral parameter equation of the pencil beam from the i-th intersection in the global coordinate system. The spiral parameter equation of the pencil beam from the i-th intersection in the global coordinate system is used to describe the theoretical motion trajectory of the pencil beam from the i-th intersection.

9. The calculation method according to claim 8, wherein: Obtaining the position of the (i+1)th intersection point according to the theoretical motion trajectory of the pencil beam from the (i)th intersection point and the boundary condition of the (i)th candidate box includes: Taking three planes in the i-th candidate box that have a probability of intersecting with the pencil beam as candidate planes; Substituting the boundary conditions of the three candidate planes into the spiral parametric equation of the pencil beam starting from the i-th intersection in the global coordinate system, a unique valid solution is obtained; Substitute the effective solution back into the spiral parameter equation of the pencil beam starting from the i-th intersection in the global coordinate system to obtain the position of the i+1-th intersection.

10. The calculation method according to claim 1, wherein: The geometric center of the phantom coincides with the geometric center of the magnetic field.

Images