Method for calculating proton radiation dose in magnetic field

Abstract

The present disclosure provides a method for calculating a proton radiation dose in a magnetic field, comprising: determining a position of an exit source and a velocity and energy of a pencil beam at the exit source; upon a proton beam being emitted from the exit source, determining, in chronological order, a plurality of candidate boxes according to the position of the exit source, the velocity and energy of the pencil beam at the exit source, and information of the magnetic field, and sequentially calculating a position of an intersection point between the pencil beam and a surface of each candidate box, and the velocity and energy of the pencil beam at each intersection point; among the candidate boxes inside a phantom, selecting a plurality of intersection points as sampling points, and connecting all the sampling points to obtain an actual motion trajectory of the pencil beam in the phantom; and according to the actual motion trajectory and a pencil beam algorithm, obtaining a proton dose delivered by the pencil beam to any voxel in the phantom, and then obtaining a proton dose delivered by proton irradiation to any voxel in the phantom.

Description

Method for calculating proton radiation dose in magnetic field

[0001] The present disclosure claims priority to the Chinese patent application No. 202410740627.4, filed on June 9, 2024, and entitled "Method for calculating proton radiation dose in magnetic field", the entire content of which is incorporated herein by reference. TECHNICAL FIELD

[0002] The present disclosure belongs to the technical field of nuclear medical imaging, and specifically relates to a method for calculating proton radiation dose in a magnetic field. BACKGROUND

[0003] Proton therapy is a new type of radiotherapy method for cancer treatment using high-energy proton beams. It has the advantages of relatively low incident dose and sharp dose falloff, and can accurately deliver dose to the target area and reduce damage to the surrounding normal tissues. However, due to factors such as uncertainty of proton range, positioning error and tissue movement, the physicist needs to increase the margin around the target area in the radiotherapy plan, which limits the high-dose gradient advantage of proton therapy.

[0004] Magnetic resonance imaging (MRI) has the advantages of real-time imaging and high soft tissue contrast, and is free of 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, so far, MRI-guided proton therapy (MRgPT) has not been realized clinically. During the dose calculation, optimization and delivery process, the dose distortion of the proton beam caused by the magnetic field effect (i.e. Lorentz force) is one of the many obstacles to realizing MRgPT.

[0005] Wolf and Bortfeld first derived the analytical solution of the proton beam deflection trajectory in a transverse uniform magnetic field. Their method is based on a small-angle approximation, so the error increases when the proton deflection angle is large; in addition, their method does not have a closed-form solution in the relativistic case. Schellhammer and Hoffmann introduced an analytical iterative method, which gives very good agreement with Monte Carlo simulation results in the energy range of 60 to 250 MeV. However, the applicability of this method is limited to the case where the magnetic field is perpendicular to the proton beam, and it is not suitable for voxelized phantoms. Padilla-Cabal proposed a numerical iterative method that stores the energy deposition values of the proton beam in water along the trajectory into a lookup table. For non-water materials, this method determines the water equivalent depth from the lookup table and scales it using a correction factor specific to the material. However, the correction factor cannot accurately describe the properties of high-density materials such as bone, leading to deviations in the dose distribution in these materials. Monte Carlo simulation provides high-precision simulation, but the calculation speed is slow, limiting its application in real-time MRgPT. Therefore, there is a need for an efficient and relatively accurate dose calculation algorithm to better achieve the calculation of radiation dose in MRgPT.

[0006] SUMMARY

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

[0008] As a first aspect of the present disclosure, a method for calculating the 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 radiate protons to the phantom, the proton beam includes a plurality of pencil beams, and the method comprises:

[0009] determining the position of the emission source and the speed and energy of the pencil beam at the emission source;

[0010] after the proton beam is emitted from the emission source, determining a plurality of candidate boxes in time sequence according to the position of the emission source, the speed and energy of the pencil beam at the emission source, and the information of the magnetic field, and sequentially calculating the position of the intersection of the pencil beam and the surface of each candidate box and the speed and energy of the pencil beam at each intersection, the candidate box corresponding to the position of the intersection located outside the phantom or the energy at the intersection being equal to zero being the last candidate box, wherein in the two adjacent candidate boxes, the candidate box set later is determined according to the position of the intersection of the pencil beam and the candidate box set earlier and the speed of the pencil beam at the intersection;

[0011] in the candidate boxes inside the phantom, selecting a plurality of intersections as sampling points, and connecting all the sampling points to obtain the actual motion trajectory of the pencil beam in the phantom;

[0012] According to the actual motion trajectory and the pencil beam algorithm, the proton dose of the pencil beam transmitted to any voxel in the phantom is obtained, and then the proton dose of the proton beam radiated to any voxel in the phantom is obtained.

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

[0014] The plurality of candidate boxes and the plurality of intersections determined in time sequence are numbered using {i|0≤i≤H}; the phantom is taken as the 0th candidate box, and the position of the emission source is taken as the 0th intersection;

[0015] The ith candidate box is determined according to the CT image data of the phantom, the position of the ith intersection, the speed and energy of the pencil beam at the ith intersection, and the information of the magnetic field, and the position of the ith+1 intersection of the pencil beam and the ith candidate box and the energy and speed of the pencil beam at the ith+1 intersection are calculated, wherein when i=H, the Hth intersection is located outside the phantom, or the speed of the pencil beam at the Hth intersection is zero, and the ith+1 intersection is no longer calculated.

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

[0017] The calculation method of the position of the ith+1 intersection of the pencil beam and the surface of the ith candidate box includes:

[0018] The theoretical motion trajectory of the pencil beam from the ith intersection is determined according to the position of the ith intersection, the speed of the pencil beam at the ith intersection, and the information of the magnetic field;

[0019] The position of the ith+1 intersection is obtained according to the theoretical motion trajectory of the pencil beam from the ith intersection and the boundary condition of the ith candidate box.

[0020] According to an embodiment of the present disclosure, the theoretical motion trajectory is a spiral trajectory, and the theoretical motion trajectory of the pencil beam from the ith intersection is determined according to the position of the ith intersection, the speed of the pencil beam at the ith intersection, and the information of the magnetic field, including:

[0021] The initial phase angle of the spiral trajectory of the pencil beam from the ith intersection is obtained according to the position of the ith intersection;

[0022] The radius and pitch of the spiral trajectory of the pencil beam from the ith intersection are determined according to the speed of the pencil beam at the ith intersection, the information of the magnetic field direction, and the energy of the pencil beam at the ith intersection;

[0023] According to the initial phase angle, radius and pitch of the spiral trajectory of the pencil beam from the i-th intersection point, a theoretical motion trajectory of the pencil beam from the i-th intersection point is obtained.

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

[0025] According to an embodiment of the present disclosure, the method for determining the energy and velocity of the pencil beam at the i+1-th intersection point comprises:

[0026] According to the CT image data of the phantom, Bragg-Kleeman parameters of the i-th candidate box in the phantom are obtained.

[0027] According to the Bragg-Kleeman parameters of the i-th candidate box, the velocity and energy of the pencil beam at the i-th intersection point, the energy and velocity of the pencil beam at the i+1-th intersection point are determined.

[0028] According to an embodiment of the present disclosure, according to the CT image data of the phantom, the Bragg-Kleeman parameters of the i-th candidate box in the phantom comprise:

[0029] The size of the voxel of the CT image is made the same as that of the i-th candidate box, wherein i>0; according to the HU value of the voxel in the CT image, the Bragg-Kleeman parameters of the i-th candidate box are extracted from a Bragg-Kleeman parameter lookup table;

[0030] According to an embodiment of the present disclosure, the magnetic field is a uniform magnetic field, and according to the initial phase angle, radius and pitch of the spiral trajectory of the pencil beam from the i-th intersection point, a theoretical motion trajectory of the pencil beam from the i-th intersection point is obtained, comprising:

[0031] According to the initial phase angle, radius and pitch of the spiral trajectory of the pencil beam from the i-th intersection point, a spiral parameter equation of the pencil beam from the i-th intersection point is established in a local coordinate system;

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

[0033] According to an embodiment of the present disclosure, 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, the position of the i+1-th intersection point is obtained, comprising:

[0034] Three planes in the i-th candidate box that have intersection probability with the pencil beam are taken as candidate planes;

[0035] The boundary conditions of the three candidate planes are brought into the spiral parameter equation of the pencil beam from the i-th intersection point in the global coordinate system to obtain a unique effective solution.

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

[0037] According to an embodiment of the present disclosure, the geometric center of the phantom coincides with the geometric center of the magnetic field. BRIEF DESCRIPTION OF DRAWINGS

[0038] FIG. 1 shows a calculation method of proton radiation dose in a magnetic field according to an embodiment of the present disclosure;

[0039] FIG. 2 shows a schematic diagram of an intersection point according to an embodiment of the present disclosure;

[0040] FIG. 3 shows a schematic diagram of a theoretical trajectory of a pencil beam in a magnetic field in a local coordinate system according to an embodiment of the present disclosure;

[0041] FIG. 4 shows a dose distribution of a proton beam incident on a prostate site at a gantry angle of 68° / 248° in a 3.0T uniform transverse magnetic field according to an embodiment of the present disclosure;

[0042] FIG. 5 shows a dose distribution of a proton beam incident on a water phantom at a gantry angle of 270° in a 3.0T uniform transverse magnetic field according to an embodiment of the present disclosure. DETAILED DESCRIPTION

[0043] In order to make the objectives, technical solutions and advantages of the present disclosure clearer, the present disclosure is further described in detail below with reference to the embodiments and the accompanying drawings. Obviously, the described embodiments are part of the embodiments of the present disclosure, rather than all the embodiments. Based on the embodiments of the present disclosure, all other embodiments obtained by those skilled in the art without creative work fall within the scope of protection of the present disclosure.

[0044] FIG. 1 shows a calculation method of proton radiation dose in a magnetic field according to an embodiment of the present disclosure. In the method, 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 radiate protons to the phantom, the proton beam includes a plurality of pencil beams, and the calculation method includes steps S1-S4.

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

[0046] In step S2, after the proton beam is emitted from the emission source, a plurality of candidate boxes are determined in time sequence according to the position of the emission source, the speed and energy of the pencil beam at the emission source, and the information of the magnetic field, and the positions of the intersection points of the pencil beam and the surface of each candidate box and the speed and energy of the pencil beam at each intersection point are calculated in sequence, and the candidate box corresponding to the position of the intersection point outside the phantom or the energy at the intersection point equal to zero is the last candidate box, wherein in the two adjacent candidate boxes, the candidate box set later is determined according to the position of the intersection point of the pencil beam and the candidate box set earlier and the speed of the pencil beam at the intersection point.

[0047] In step S3, in the candidate boxes inside the phantom, a plurality of 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.

[0048] In step S4, the proton dose of the pencil beam transferred to any voxel in the phantom is obtained according to the actual motion trajectory and the pencil beam algorithm, and then the proton dose of the pencil beam radiated to any voxel in the phantom is obtained.

[0049] According to the embodiments of the present disclosure, the actual motion trajectory of the pencil beam in the phantom is determined by setting a plurality of candidate boxes and calculating the intersection points of the pencil beam and the surface of each candidate box, and the proton dose of the pencil beam transferred to any voxel in the phantom is obtained in combination with the actual motion trajectory and the pencil beam algorithm, so that the determination of the proton radiation dose in the magnetic field is realized. The method provided by the embodiments of the present disclosure has high precision. The method of the present disclosure can be used to realize the magnetic resonance guided proton therapy for patients.

[0050] According to the embodiments of the present disclosure, in step S2, after the proton beam is emitted from the emission source, a plurality of candidate boxes are determined in time sequence according to the position of the emission source, the speed and energy of the pencil beam at the emission source, and the information of the magnetic field, and the positions of the intersection points of the pencil beam and the surface of each candidate box and the speed and energy of the pencil beam at each intersection point are calculated in sequence, and the candidate box corresponding to the position of the intersection point outside the phantom or the energy at the intersection point equal to zero is the last candidate box, wherein in the two adjacent candidate boxes, the candidate box set later is determined according to the position of the intersection point of the pencil beam and the candidate box set earlier and the speed of the pencil beam at the intersection point.

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

[0052] In step S22, the i th candidate box is determined according to the CT image data of the phantom, the position of the i th intersection point, the speed and energy of the pencil beam at the i th intersection point, and the information of the magnetic field, and the position of the i+1 th intersection point of the pencil beam and the i th candidate box and the energy and speed of the pencil beam at the i+1 th intersection point are calculated, wherein when i=H, the H th intersection point is located outside the phantom, or the speed of the pencil beam at the H th intersection point is zero, and the i+1 th intersection point is no longer calculated.

[0053] FIG. 2 shows a schematic diagram of intersection points according to an embodiment of the present disclosure.

[0054] As shown in FIG. 2, P i+1 The arrow direction represents the velocity direction of the pencil beam at each intersection point.

[0055] According to an embodiment of the present disclosure, in step S1, the exit source can be located in the magnetic field or outside the magnetic field, which is related to the size of the magnetic field. The exit source in step S1 is not a physical point source in the head of the treatment machine. The proton beam is generated at the physical point source, passes through various components in the head of the proton therapy machine, and is scanned and diffused into a wider two-dimensional plane, which is the exit source. The exit source is located at the trailing edge of the range shifter, and the exit source can also be called a virtual source.

[0056] According to an embodiment of the present disclosure, the velocity of the pencil beam at the exit source includes the size and direction, the velocity and energy of the pencil beam at the position of the exit source, and the position of the exit source can be obtained according to the radiotherapy plan file. The information provided by the radiotherapy plan file includes the number of irradiation points (i.e., the number of pencil beams covering the phantom), the gantry angle, the head position, the isocenter, the source axis distance, the range shifter information of the treatment machine in the beam eye coordinate system, the treatment bed angle for placing the phantom, and the irradiation time (in units of MU) of each control point (pencil beam), the nominal energy of the pencil beam. The initial phase space information can be obtained according to the radiotherapy plan 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 position of the exit source), the divergence angle, the beam spot size, and the correlation coefficient. The energy of the pencil beam at the position of the exit source can be obtained by using the beam energy model.

[0057] According to an embodiment of the present disclosure, the angle between the velocity direction of the pencil beam and the direction of the magnetic field is θ, and θ ∈ [0, π]. In step S22, the calculation method of the position of the i+1 intersection point where the pencil beam intersects the surface of the i candidate box includes steps S221-S222.

[0058] In step S221, the theoretical motion trajectory of the pencil beam from the i intersection point is determined according to the position of the i intersection point, the velocity of the pencil beam at the i intersection point, and the information of the magnetic field.

[0059] In step S222, the position of the i+1 intersection point is obtained according to the theoretical motion trajectory of the pencil beam from the i intersection point and the boundary condition of the i candidate box.

[0060] According to an embodiment of the present disclosure, the above-mentioned theoretical motion trajectory is a spiral trajectory, and in step S222,

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

[0062] In step S2211, an initial phase angle of a spiral trajectory of the pencil beam from the i-th intersection is obtained according to the position of the i-th intersection.

[0063] In step S2212, a radius and a pitch of the spiral trajectory of the pencil beam from the i-th intersection are determined according to the speed of the pencil beam at the i-th intersection, the information of the direction of the magnetic field, and the energy of the pencil beam at the i-th intersection.

[0064] In step S2213, the theoretical motion trajectory of the pencil beam from the i-th intersection is obtained according to the initial phase angle, the radius, and the pitch of the spiral trajectory of the pencil beam from the i-th intersection.

[0065] According to an embodiment of the present disclosure, the magnetic field described above can be a uniform magnetic field or a non-uniform magnetic field, for example, when the magnetic field is a uniform magnetic field, the magnetic field strength can be 3.0T, the direction of the magnetic field can be along the negative direction of the Z axis of the IEC coordinate system, and the shape can be a cuboid (the length, width, and height can all be 50cm). In step S2213, the theoretical motion trajectory of the pencil beam from the i-th intersection is obtained according to the initial phase angle, the radius, and the pitch of the spiral trajectory of the pencil beam from the i-th intersection, including steps S22131-S22132.

[0066] In step S22131, a spiral parameter equation of the pencil beam from the i-th intersection is established in the local coordinate system by using the initial phase angle, the radius, and the pitch of the spiral trajectory of the pencil beam from the i-th intersection.

[0067] According to an embodiment of the present disclosure, without considering the energy loss for the moment, the spiral parameter equation P lcs (φ) in the local coordinate system (LCS) is as follows. lcs (φ), y lcs (φ), z lcs (φ)) is as follows.

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

[0069] In step S22132, the spiral parameter equation of the pencil beam from the i-th intersection is converted to the global coordinate system to obtain a spiral parameter equation of the pencil beam from the i-th intersection in the global coordinate system, which is used to describe the theoretical motion trajectory of the pencil beam from the i-th intersection.

[0070] An orthogonal matrix M composed of LGS base vectors i , P lcs is converted to GCS to obtain the running trajectory equation P(φ) of the proton in GCS. According to an embodiment of the present disclosure, the parametric equation (formula (1)) of the spiral trajectory can only be expressed in a local coordinate system. The Z-axis of the local coordinate system is anti-parallel to the magnetic field direction, and the orthogonal base matrix of the local coordinate system is denoted as M i

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

[0072] In step S223, the Bragg-Kleeman parameters of the i-th candidate box in the phantom are obtained according to the CT image data of the phantom.

[0073] In step S224, the energy and speed of the pencil beam at the i+1-th intersection point 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 point.

[0074] According to an embodiment of the present disclosure, in 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 includes steps S2231-S2235.

[0075] In step S2231, the voxels of the CT image are made the same size as the i-th candidate box, where i>0.

[0076] In step S2232, the Bragg-Kleeman parameters of the i-th candidate box are extracted from the Bragg-Kleeman parameter lookup table according to the HU values of the voxels in the CT image.

[0077] According to an embodiment of the present disclosure, the curvature of the spiral motion trajectory of the pencil beam in the magnetic field is related to the energy thereof, in order to maximize the accuracy of the deflection prediction, the embodiment of the present disclosure uses a contrast table of HU value and Bragg-Kleeman parameters (a, p), material mass density, so as to consider the influence of the energy change on the relative stopping power (RSP) of the proton in the calculation process; after the HU value in the range of-1000-2995 is converted into 25 materials according to the parameter conversion file, then the integral depth dose curve (IDD) of the pencil beam with the energy range of 10-250 MeV and the energy interval of 10 MeV in the 25 materials is simulated by using the Monte Carlo program, so as to obtain the range R under each nominal energy E, and the fitting of E and R is performed, that is, the Bragg-Kleeman parameters of each material are obtained. The embodiment of the present disclosure can also select to use the HU-RSP calibration curve to calculate the beam trajectory.

[0078] In step S224, the energy and the speed of the pencil beam at the i+1 intersection point are determined according to the Bragg-Kleeman parameters of the i candidate box, the speed and the energy of the pencil beam at the i intersection point, and can be expressed as follows. According to an embodiment of the present disclosure, the medium will consume the energy of the proton beam, so the residual energy E of the pencil beam at the i+1 intersection point is corrected by using the Bragg-Kleeman rule i+1 :

[0079] According to an embodiment of the present disclosure, a i , p i is the Bragg-Kleeman parameter of the i candidate box, p i is the mass density of the i candidate box, and l arc is the distance walked by the pencil beam in the i candidate box.

[0080] According to an embodiment of the present disclosure, due to the influence of the Lorentz force, the speed direction of the proton beam perpendicular to the magnetic field direction will be deflected, so it is necessary to correct the speed direction by using the rotation matrix R z rotating around the Z axis. If the proton is in the air outside the medium, only the direction is corrected, and if it is in the voxel, the energy is also corrected. The speed direction correction formula is as follows:

[0081] wherein, is the unit vector of the speed of the pencil beam at the i+1 intersection point, is the component of the unit vector of the speed of the pencil beam at the i intersection point in parallel to the magnetic field direction, is the component of the unit vector of the speed of the pencil beam at the i intersection point perpendicular to the magnetic field direction, is the transpose matrix of M i , and Rz (φ i+1 ) T for R z (φ i+1 ) is the transpose matrix of φ i+1 , φ

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

[0083] When the proton source is located outside the magnetic field, the theoretical trajectory of the pencil beam before reaching the phantom is determined by the position where the proton source is located, the geometric boundary of the magnetic field, the direction and strength of the magnetic field, the energy and direction of the pencil beam. This is because the pencil beam moves along a straight line before entering the magnetic field, and starts to move along a curve after entering the magnetic field due to the Lorentz force. Therefore, when the proton source is located outside the magnetic field, the trajectory of the pencil beam before entering the magnetic field cannot be determined using the spiral trajectory.

[0084] According to an embodiment of the present disclosure, in step S222, the position of the i+1 intersection point is obtained according to the theoretical trajectory of the pencil beam from the i intersection point and the boundary condition of the i candidate box, comprising steps S2221-S2223.

[0085] In step S2221, three planes in the i candidate box that have intersection probability with the pencil beam are taken as candidate planes.

[0086] In step S2222, the boundary conditions of the three candidate planes are brought into the spiral parameter equation of the pencil beam from the i intersection point in the global coordinate system to obtain a unique valid solution.

[0087] In step S2223, the valid solution is substituted back into the spiral parameter equation of the pencil beam from the i intersection point in the global coordinate system to obtain the position of the i+1 intersection point.

[0088] According to an embodiment of the present disclosure, the geometric body representing the phantom is taken as the 0 candidate box, and the three surfaces opposite to the momentum direction of the protons in the pencil beam are taken as the candidate planes; the voxel i where the pencil beam is located is taken as the i candidate box, and the three surfaces opposite to the momentum direction of the protons in the pencil beam are taken as the candidate planes.

[0089] According to an embodiment of the present disclosure, for the 0 candidate box, the solution interval is set to [0, 0.3π], and for the i candidate box (i>0), the solution interval is [0, the diagonal length of the i candidate box / a i ]; by judging whether the intersection point P i+1 is located on the surface of the candidate box, a unique valid solution φ i+1 is obtained from the three solutions.

[0090] According to an embodiment of the present disclosure, in step S3, the line between two adjacent sampling points in the n sampling points is a straight line, 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 interval of the sampling points can be adjusted according to the requirements of calculation accuracy and efficiency.

[0091] According to an embodiment of the present disclosure, after step S3, the above method further comprises: calculating the physical path length l and the radial distance r of each voxel in the phantom. The radial distance is the perpendicular line of the actual motion trajectory passing through the center of a voxel, and the distance from the center of the voxel to the foot of the perpendicular line. The physical path length is the distance from the first sampling point to the foot of the perpendicular line along the actual motion trajectory of the pen-shaped beam. For an actual motion trajectory, each voxel has a radial distance and a physical path length.

[0092] According to an embodiment of the present disclosure, in step S4, the proton motion trajectory is combined with the pen-shaped beam algorithm, and the dose d(r, l, l w ) delivered by a pen-shaped beam to a voxel located at the point of interest (r, l) can be calculated by multiplying the integral depth dose (IDD) by the double Gaussian kernel function (K), using the formula as shown below:

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

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

[0095] wherein 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 reaction, σ mcs is the variance of the multiple Coulomb scattering Gaussian kernel, and σ nuc is the variance of the nuclear reaction Gaussian kernel. Wherein σ mcs and σ nuc are obtained by the divergence angle, the beam spot size and the correlation coefficient in the initial phase space information.

[0096] wherein IDD(l w) are extracted from the IDD database according to the equivalent water depth. The IDD database has a depth range of 0-40 cm, a resolution of 0.1 mm, an energy range of 10.1-250.0 MeV, and an energy resolution of 0.1 MeV. The IDD database is established by a Monte Carlo simulation method.

[0097] According to an embodiment of the present disclosure, the dose of all pencil beams (PBs) is superimposed to obtain a dose distribution D(x, y, z) of a voxel located at (x, y, z) in a magnetic field.

[0098] wherein w PB is the weight of each pencil beam.

[0099] According to an embodiment of the present disclosure, a Gaussian energy spectrum is obtained according to IDD measurement data of pencil beams emitted by an actual treatment machine at different nominal energies, from which the energy spread and the center energy of the beam can be extracted. According to air ionization (IAF) measurement data of the treatment machine, the initial beam spot size, angular divergence, and correlation of the pencil beam are obtained. The pencil beam is detected by a detector, and according to the dose deposition in the detector, the number of emitted particles per 1 MU of the pencil beam at different nominal energies is determined. The number of emitted particles per 1 MU of the pencil beam is used to determine the weight of the pencil beam.

[0100] According to an embodiment of the present disclosure, the determination of the proton radiation dose in a magnetic field can also be based on an MRI image, but since the MRI signal intensity depends on the proton density and its tissue relaxation characteristics, it cannot be directly used for dose calculation, and the MRI needs to be converted into a CT HU image, i.e., a "pseudo-CT (sCT)" image. A mapping function of the MRI voxel intensity and the HU value is constructed by using a large amount of clinical data to realize the conversion of the MRI intensity value and the HU value, and a radiotherapy plan of the ion beam is made based on the sCT.

[0101] FIG. 4 is a dose distribution of a proton beam incident on a prostate site at a gantry angle of 68° / 248° under a 3.0T uniform transverse magnetic field.

[0102] Part (a) of FIG. 4 is a result determined by using an embodiment of the present disclosure, part (b) is a simulation result by a Monte Carlo method, and part (c) is a 2mm / 2% Gamma pass rate map of the result of the present disclosure. As can be seen from FIG. 4, the consistency with the Monte Carlo simulation result shows the accuracy of the present disclosure, and the 2mm / 2% Gamma pass rate at a 10% threshold is 99.20%.

[0103] FIG. 5 is a dose distribution of a proton beam incident on a water phantom at a gantry angle of 270° under a 3.0T uniform transverse magnetic field.

[0104] The result determined by the embodiment of the present disclosure in (a) of FIG. 5, the simulation result of the Monte Carlo method in (b), and the 2mm / 2% Gamma passing rate diagram in (c) of the result of the present disclosure. As can be seen from FIG. 5, the consistency with the Monte Carlo simulation result illustrates the accuracy of the present disclosure, and the 2mm / 2% Gamma passing rate at a 10% threshold is 98.97%.

[0105] In actual application, for example, when performing proton therapy on a patient, an organ sketch file, a CT image of the patient, and a radiotherapy plan file need to be obtained.

[0106] The embodiment of the present disclosure uses a local coordinate system and a global coordinate system, derives a theoretical trajectory equation of a pencil beam in a magnetic field based on a spiral parameter equation and a Lorentz equation under relativity, substitutes the boundary conditions of the candidate box surface, numerically solves the motion trajectory of the pencil beam in the phantom under the magnetic field, and finally combines the pencil beam algorithm to calculate the dose distribution of proton radiation under the magnetic field, thereby solving three main problems in the field. First, the applicability of the existing analytical algorithm for ray tracing of charged particles in a magnetic field is limited to the case where the magnetic field is perpendicular to the proton beam, and is not applicable to voxelized phantoms. Second, the existing numerical iterative algorithm for charged particles in a magnetic field cannot accurately describe the properties of high-density materials, resulting in large errors in the dose distribution in these materials. Third, the existing traditional Monte Carlo program is slow in calculation and cannot meet the speed requirements of real-time dose calculation in actual MRgRT. The present disclosure solves a series of problems described above and can be used for the calculation of proton radiation dose in a magnetic field in a treatment planning system, shortens the calculation time while ensuring the calculation accuracy.

[0107] The method provided by the embodiment of the present disclosure 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.

[0108] The method provided by the embodiment of the present disclosure only generates a single-variable nonlinear function when calculating the actual motion trajectory of the pencil beam, and the solving range can be set according to the size of the candidate box, so that the number of iterations required for numerical solving is greatly reduced, and the calculation speed is improved.

[0109] The method provided by the embodiment of the present disclosure approximates the proton trajectory as a polygonal line connected by sampling points by sampling the intersection points, to calculate the physical path length l and the radial distance r of each voxel in the medium, and successfully combines with the pencil beam algorithm.

[0110] The above specific embodiments further illustrate the purpose, technical solutions and advantages of the present disclosure. It should be understood that the above are only specific embodiments of the present disclosure and are not used to limit the present disclosure. Any modification, equivalent replacement, improvement, etc. made within the spirit and principles of the present disclosure shall be included in the protection scope of the present disclosure.

Claims

1. A method of calculating the dose of proton radiation in a magnetic field, wherein, Proton beams are emitted from an emission source and enter a magnetic field, the proton beams move in the magnetic field and enter a phantom in the magnetic field to radiate protons to the phantom, the proton beams include a plurality of pencil beams, and the method comprises: determining the position of the emission source and the speed and energy of the pencil beams at the emission source; after the proton beams are emitted from the emission source, a plurality of candidate boxes are determined in time sequence according to the position of the emission source, the speed and energy of the pencil beams at the emission source, and the information of the magnetic field, and the positions of the intersection points of the pencil beams with the surface of each candidate box and the speed and energy of the pencil beams at each intersection point are calculated in turn, the last candidate box is the candidate box corresponding to the position of the intersection point outside the phantom or the energy at the intersection point being equal to zero, wherein in the two adjacent candidate boxes, the candidate box arranged later is determined according to the position of the intersection point of the pencil beams with the candidate box arranged earlier and the speed of the pencil beams at the intersection point; in the candidate boxes inside the phantom, a plurality of intersection points are selected as sampling points, and all the sampling points are connected to obtain the actual motion trajectory of the pencil beams in the phantom; according to the actual motion trajectory and the pencil beam algorithm, the proton dose delivered by the pencil beams to any voxel in the phantom is obtained, and then the proton dose radiated by the pencil beams to any voxel in the phantom is obtained.

2. The calculation method according to claim 1, wherein after the proton beams are emitted from the emission source, a plurality of candidate boxes are determined in time sequence according to the position of the emission source, the speed and energy of the pencil beams at the emission source, and the information of the magnetic field, and the positions of the intersection points of the pencil beams with the surface of each candidate box and the speed and energy of the pencil beams at each intersection point are calculated in turn, which comprises: using {i|0≤i≤H} to number the plurality of candidate boxes and the plurality of intersection points determined in time sequence; taking the phantom as the 0th candidate box and the emission source position as the 0th intersection point; determining the ith candidate box according to the CT image data of the phantom, the position of the ith intersection point, the speed and energy of the pencil beams at the ith intersection point, and the information of the magnetic field, and calculating the position of the ith+1 intersection point of the pencil beams with the ith candidate box and the energy and speed of the pencil beams at the ith+1 intersection point, wherein when i=H, the Hth intersection point is located outside the phantom, or the speed of the pencil beams at the Hth intersection point is zero, and the ith+1 intersection point is no longer calculated.

3. The computational method of claim 2, wherein, the angle between the speed direction of the pencil beam and the direction of the magnetic field is θ, and θ∈[0, π]; the calculation method of the position of the ith+1 intersection point of the pencil beams intersecting with the surface of the ith candidate box comprises: determining the theoretical motion trajectory of the pencil beams from the ith intersection point according to the position of the ith intersection point, the speed of the pencil beams at the ith intersection point, and the information of the magnetic field; obtaining the position of the ith+1 intersection point according to the theoretical motion trajectory of the pencil beams from the ith intersection point and the boundary condition of the ith candidate box.

4. The computational method of claim 3, wherein, the theoretical motion trajectory is a spiral trajectory, determining a theoretical trajectory of the pencil beam from the ith intersection point according to the position of the ith intersection point, the velocity of the pencil beam at the ith intersection point, and information of the magnetic field, comprising: obtaining an initial phase angle of a helical trajectory of the pencil beam from the ith intersection point according to the position of the ith intersection point; determining a radius and a pitch of the helical trajectory of the pencil beam from the ith intersection point according to the velocity of the pencil beam at the ith intersection point, information of the direction of the magnetic field, and the energy of the pencil beam at the ith intersection point; obtaining the theoretical trajectory of the pencil beam from the ith intersection point according to the initial phase angle, the radius, and the pitch of the helical trajectory of the pencil beam from the ith intersection point.

5. The calculation method of claim 2, further comprising, at each intersection point in the magnetic field, correcting the velocity of the pencil beam by using a coordinate base matrix.

6. The computational method of 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 ith candidate box in the phantom according to CT image data of the phantom; determining the energy and velocity of the pencil beam at the (i+1)th intersection point according to the Bragg-Kleeman parameters of the ith candidate box, the velocity, and the energy of the pencil beam at the ith intersection point.

7. The computational method of claim 6, wherein, The method for obtaining Bragg-Kleeman parameters of the ith candidate box in the phantom according to CT image data of the phantom comprises: making the size of a voxel of the CT image the same as that of the ith candidate box, wherein i>0; extracting the Bragg-Kleeman parameters of the ith candidate box in a Bragg-Kleeman parameter lookup table according to the HU value of the voxel in the CT image.

8. The computational method of claim 3, wherein, The magnetic field is a uniform magnetic field, and the method for obtaining the theoretical trajectory of the pencil beam from the ith intersection point according to the initial phase angle, the radius, and the pitch of the helical trajectory of the pencil beam from the ith intersection point comprises: establishing a helical parameter equation of the pencil beam from the ith intersection point in a local coordinate system by using the initial phase angle, the radius, and the pitch of the helical trajectory of the pencil beam from the ith intersection point; converting the helical parameter equation of the pencil beam from the ith intersection point in the local coordinate system to a helical parameter equation of the pencil beam from the ith intersection point in a global coordinate system, wherein the helical parameter equation of the pencil beam from the ith intersection point in the global coordinate system is applicable to describing the theoretical trajectory of the pencil beam from the ith intersection point.

9. The computational method of claim 3, wherein, The method for obtaining the position of the (i+1)th intersection point according to the theoretical trajectory of the pencil beam from the ith intersection point and the boundary condition of the ith candidate box comprises: taking three planes in the ith candidate box that have intersection probability with the pencil beam as candidate planes; bringing the boundary conditions of the three candidate planes into the helical parameter equation of the pencil beam from the ith intersection point in the global coordinate system to obtain a unique valid solution; substituting the valid solution back into the helical parameter equation of the pencil beam from the ith intersection point in the global coordinate system to obtain the position of the (i+1)th intersection point.

10. The computational method of claim 1, wherein, The geometric center of the phantom coincides with the geometric center of the magnetic field.

Images