Graphical processor-based proton dose distribution determination method and apparatus

Abstract

The application provides a kind of based on graphics processor's proton dose distribution determination method and device, can be applied to medical radiation dose calculation technical field.The method comprises: calculating the respective pen beam trajectory of multiple pen beams;Obtain the respective dose kernel parameter of pen beam at multiple intersection points;According to field dose cutoff radius, determine the voxel of interest from the tissue voxel in the tissue voxel grid, wherein the field dose cutoff radius represents the upper limit value of the action distance of all pen beams under one field that significantly contributes to dose distribution in the transverse direction;Based on the respective dose kernel parameter of multiple intersection points, call multiple threads and execute multiple dose calculation tasks of voxel of interest in parallel, obtain the voxel deposition dose of voxel of interest;According to the respective voxel deposition dose of multiple voxel of interest, determine the dose distribution result related to specified tissue.

Description

Technical Field

[0001] This invention relates to the field of medical radiation dose calculation technology, and specifically to a method and apparatus for determining proton dose distribution based on a graphics processor. Background Technology

[0002] MRI-guided proton therapy (MRPT) is an advanced radiotherapy technique that integrates magnetic resonance imaging (MRI) with proton therapy. MRPT leverages the real-time imaging capabilities, high soft-group contrast, and lack of ionizing radiation provided by MRI to offer precise image guidance for proton dose delivery, thereby delivering the dose more accurately to the patient's tumor tissue and maximizing the protection of healthy tissue. As MRPT evolves towards adaptive radiation therapy (ART), it places increasingly stringent demands on the efficiency of dose calculation. However, current GPU-accelerated Monte Carlo algorithms still struggle to meet the computational efficiency requirements of real-time applications, while related analytical algorithms have limited accuracy in dose calculation within complex, non-uniform MRI magnetic fields. Therefore, MRPT currently lacks a dose calculation method that simultaneously achieves high efficiency and high accuracy. Summary of the Invention

[0003] In view of the above problems, the present invention provides a method and apparatus for determining proton dose distribution based on a graphics processor.

[0004] According to a first aspect of the present invention, a proton dose distribution determination method based on a graphics processor is provided, comprising: calculating the pencil beam trajectories of multiple pencil beams, wherein the intersection points in the pencil beam trajectories represent the intersection of the pencil beams with a medium voxel grid representing the transmission medium, or the intersection of the pencil beams with a tissue voxel grid representing a specified tissue; obtaining the dose kernel parameters of the pencil beams at the multiple intersection points; determining the voxel of interest from the tissue voxels in the tissue voxel grid according to the field dose cutoff radius, wherein the field dose cutoff radius represents the upper limit of the action distance of all pencil beams in a field that makes a significant contribution to the dose distribution in the lateral direction; based on the dose kernel parameters of the multiple intersection points, calling multiple threads to execute the dose calculation tasks of the multiple voxels of interest in parallel to obtain the voxel deposition dose of the voxel of interest; and determining the dose distribution result related to the specified tissue according to the voxel deposition dose of the multiple voxels of interest.

[0005] A second aspect of the present invention provides a proton dose distribution determination device based on a graphics processor, comprising: a trajectory calculation module for calculating the trajectories of multiple pencil beams, wherein intersection points in the pencil beam trajectories represent intersections between the pencil beams and a medium voxel grid representing the transmission medium, or intersections between the pencil beams and a tissue voxel grid representing a specified tissue; an acquisition module for acquiring dose kernel parameters of the pencil beams at the multiple intersection points; a first determination module for determining voxels of interest from tissue voxels in the tissue voxel grid based on the field dose cutoff radius, wherein the field dose cutoff radius represents the upper limit of the action distance of all pencil beams in a field that significantly contribute to the dose distribution in the lateral direction; a calling module for calling multiple threads to execute dose calculation tasks of multiple voxels of interest in parallel based on the dose kernel parameters of the multiple intersection points, to obtain the voxel deposition dose of the voxels of interest; and a second determination module for determining the dose distribution result related to a specified tissue based on the voxel deposition dose of the multiple voxels of interest.

[0006] A third aspect of the present invention provides an electronic device comprising: one or more processors; and a memory for storing one or more programs, wherein when the one or more programs are executed by the one or more processors, the one or more processors perform the method described above.

[0007] According to the proton dose distribution determination method and apparatus based on a graphics processor provided by the present invention, since the field dose cutoff radius represents the upper limit of the action distance of all pencil beams in a field that significantly contribute to the dose distribution in the lateral direction, the voxels of interest are screened in the tissue voxel grid based on the field dose cutoff radius. Based on this limitation of the region of interest, the dose superposition calculation is strictly constrained to the area directly irradiated by the pencil beam and its surrounding boundary extension area smaller than the field dose cutoff radius, thereby significantly reducing redundant calculations and significantly improving the dose calculation efficiency. In addition, a voxel-level parallel strategy is adopted to realize the parallel processing of dose calculation tasks for multiple voxels of interest. By making full use of the large-scale parallel architecture of the graphics processor, the dose calculation task is decomposed into multiple independently executable dose calculation subtasks, and a large number of threads are scheduled to calculate the voxel deposition dose of the voxels of interest concurrently, further improving the dose calculation efficiency. Attached Figure Description

[0008] The above and other objects, features and advantages of the present invention will become clearer from the following description of embodiments of the invention with reference to the accompanying drawings.

[0009] Figure 1 A schematic diagram of the target object placement in the MRPT system according to an embodiment of the present invention, using parallel magnetic field configuration and vertical magnetic field configuration respectively.

[0010] Figure 2A flowchart of a proton dose distribution determination method based on a graphics processor according to an embodiment of the present invention is shown.

[0011] Figure 3 A schematic diagram illustrating a two-dimensional magnetic field distribution according to an embodiment of the present invention is shown.

[0012] Figure 4 A flowchart of a proton dose distribution determination method based on a graphics processor according to another embodiment of the present invention is shown.

[0013] Figure 5 A schematic diagram showing the dose distribution results of a simulated proton incident on a water model according to an embodiment of the present invention is illustrated.

[0014] Figure 6 A schematic diagram showing the dose distribution results of simulated proton incident on a liver tumor according to an embodiment of the present invention is illustrated.

[0015] Figure 7 A structural block diagram of a proton dose distribution determination device based on a graphics processor according to an embodiment of the present invention is shown.

[0016] Figure 8 A block diagram of an electronic device suitable for implementing a graphics processor-based proton dose distribution determination method according to an embodiment of the present invention is shown. Detailed Implementation

[0017] Hereinafter, embodiments of the present invention will be described with reference to the accompanying drawings. However, it should be understood that these descriptions are exemplary only and are not intended to limit the scope of the invention. In the following detailed description, numerous specific details are set forth to provide a thorough understanding of the embodiments of the invention for ease of explanation. However, it will be apparent that one or more embodiments may be practiced without these specific details. Furthermore, descriptions of well-known structures and techniques are omitted in the following description to avoid unnecessarily obscuring the concept of the invention.

[0018] The terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the invention. The terms “comprising,” “including,” etc., as used herein indicate the presence of features, steps, operations, and / or components, but do not exclude the presence or addition of one or more other features, steps, operations, or components.

[0019] All terms used herein (including technical and scientific terms) have the meanings commonly understood by those skilled in the art, unless otherwise defined. It should be noted that the terms used herein are to be interpreted in a manner consistent with the context of this specification, and not in an idealized or overly rigid way.

[0020] When using expressions such as "at least one of A, B and C", they should generally be interpreted in accordance with the meaning that is commonly understood by those skilled in the art (e.g., "a system having at least one of A, B and C" should include, but is not limited to, a system having A alone, a system having B alone, a system having C alone, a system having A and B, a system having A and C, a system having B and C, and / or a system having A, B and C, etc.).

[0021] In the process of developing this invention, it was discovered that MRPT has gradually become a research hotspot in the field of radiotherapy in recent years. This technology combines the advantages of both MRI and proton therapy: on the one hand, it leverages the advantages of MRI in terms of high soft tissue contrast, non-ionizing radiation, and real-time imaging; on the other hand, it leverages the dosimetric advantages of proton therapy. However, the realization of MRPT still faces challenges in the development of dose calculation engines.

[0022] MRPT dose calculation needs to balance accuracy and efficiency in realistic MRI magnetic fields. Among related methods, Monte Carlo simulation has high accuracy but is time-consuming; analytical methods are fast but lack accuracy in complex magnetic fields. The few GPU-accelerated MRPT dose calculation methods that have emerged in recent years are still limited to uniform or simplified magnetic field models, generally ignoring the spatial variation characteristics of real three-dimensional magnetic fields, thus their clinical applicability is severely limited.

[0023] Achieving accurate and efficient MRPT dose calculation faces three major challenges: First, the complexity of the MRI magnetic field. The non-uniform main magnetic field of MRI causes complex proton deflection, significantly increasing the difficulty of dose calculation. Second, the diversity of clinical configurations. Clinical MRPT systems have two main magnetic field configurations: parallel and perpendicular. Current methods have not yet been able to simultaneously accommodate both within a unified algorithmic framework. Third, the computational architecture and efficiency of current methods are bottlenecks. Current dose calculations mainly rely on serial execution by a central processing unit (CPU), and the process involves a large number of repetitive and redundant operations, resulting in a computational speed far from meeting the second-level requirements of APT.

[0024] In view of this, embodiments of the present invention provide a method and apparatus for determining proton dose distribution based on a graphics processor. The method includes: calculating the trajectory of each of multiple pencil beams; obtaining the dose kernel parameters of each pencil beam at multiple intersection points; determining the voxel of interest from tissue voxels in a tissue voxel grid based on the field dose cutoff radius, wherein the field dose cutoff radius characterizes the upper limit of the action distance of all pencil beams in a field that significantly contribute to the dose distribution in the lateral direction; based on the dose kernel parameters of each of the multiple intersection points, calling multiple threads to execute the dose calculation task of multiple voxels of interest in parallel to obtain the voxel deposition dose of the voxel of interest; and determining the dose distribution result related to a specified tissue based on the voxel deposition dose of each of the multiple voxels of interest.

[0025] In the technical solution of this invention, the user information (including but not limited to user personal information, user image information, user device information, such as location information) and data (including but not limited to data used for analysis, stored data, displayed data, etc.) involved are all information and data authorized by the user or fully authorized by all parties. Furthermore, the collection, storage, use, processing, transmission, provision, disclosure, and application of related data all comply with relevant laws, regulations, and standards, take necessary confidentiality measures, do not violate public order and good morals, and provide corresponding operation entry points for users to choose to authorize or refuse.

[0026] Figure 1 A schematic diagram of the target object placement in the MRPT system according to an embodiment of the present invention, using parallel magnetic field configuration and vertical magnetic field configuration respectively.

[0027] like Figure 1 As shown, the magnets in the split-type MRI device adopt a separate design. A physical gap exists between the first superconducting magnet unit 101 and the second superconducting magnet unit 102, providing open space for the target object region 103 and the beam path, thus distinguishing it from traditional closed-aperture structures. The thin cylinder 104 in the figure represents the treatment head, from which the pencil beam exits and enters the target object. Figure 1 a is a schematic diagram of the parallel magnetic field configuration: the target object is placed between the first superconducting magnet unit 101 and the second superconducting magnet unit 102, and the beam enters the target object through the magnet hole; Figure 1 Figure b shows a schematic diagram of a vertical magnetic field configuration: the target object is placed along the long axis of the magnet aperture, and the beam enters the target object through the gap between the magnet poles. In both magnetic field configurations, the isocenter is precisely aligned with the geometric center of the magnetic field. The "isocenter" refers to the spatial reference point formed by the intersection of the gantry rotation axis, collimator rotation axis, and treatment bed rotation axis in a radiotherapy system; this point is typically set at the geometric center of the tumor region.

[0028] To achieve multi-field irradiation, the following methods can be used: In a parallel magnetic field configuration, the MRI and proton therapy gantry are relatively fixed, and the beam injection direction is changed by synchronously rotating the magnet-gantry integrated system; in a vertical magnetic field configuration, the MRI remains fixed, and different angle irradiation is achieved by independently rotating the treatment gantry. These hardware coupling and motion methods represent two feasible technical paths for achieving multi-field irradiation in current MRPT prototypes.

[0029] It should be noted that the sequence numbers of the operations in the following methods are for descriptive purposes only and should not be considered as indicating the execution order of the operations. Unless explicitly stated otherwise, the method does not need to be executed in the exact order shown.

[0030] Figure 2 A flowchart of a proton dose distribution determination method based on a graphics processor according to an embodiment of the present invention is shown.

[0031] like Figure 2 As shown, the proton dose distribution determination method based on a graphics processor includes operations S210 to S250.

[0032] In operation S210, the individual pencil beam trajectories of multiple pencil beams are calculated.

[0033] In operation S220, the dose kernel parameters of the pencil beam at multiple intersection points are obtained.

[0034] In operation S230, the voxel of interest is determined from the tissue voxels in the tissue voxel grid based on the field dose cutoff radius.

[0035] In operation S240, based on the dose kernel parameters of each of the multiple intersection points, multiple threads are called to execute the dose calculation tasks of multiple voxels of interest in parallel, so as to obtain the voxel deposition dose of the voxels of interest.

[0036] In operation S250, dose distribution results related to a specified tissue are determined based on the voxel deposition doses of each of the multiple voxels of interest.

[0037] A pencil beam is a proton beam used in pencil beam scanning technology. The beam spot diameter is on the order of millimeters, the energy distribution is narrow and diffused, and it consists of a large number of protons. The pencil beam can be precisely guided to a predetermined target point in three-dimensional space via a deflecting magnet system in the treatment head.

[0038] The GPU-based proton dose distribution determination method is implemented using multiple CUDA (Compute Unified Device Architecture) kernel functions written in a programming language. These kernel functions are code that executes in parallel on the GPU and can run concurrently on multiple threads. First, multiple pencil beams are initialized based on the proton beam phase space source, each defined by its exit parameters in phase space (including position, direction, and energy). Then, the first CUDA kernel function is called, driving multiple threads to perform pencil beam trajectory calculations in parallel. Each thread independently handles the full ray tracing of one pencil beam.

[0039] In pencil-beam ray tracing, the intersection points of the pencil beam trajectory and the axis-aligned bounding box (AABB) mesh are calculated. This mesh is a regular mesh that discretizes the three-dimensional space, where each AABB is a cube with sides parallel to the coordinate axes. Within the area defined by the computed tomography (CT) scan, one AABB corresponds to one CT voxel, and is assigned corresponding tissue material properties based on its CT value; this part is the tissue voxel mesh. Outside this area, the algorithm automatically generates a mesh with a resolution of 1×1×1 mm. 3 A regular grid is used to cover the entire magnetic field space. This part of the grid is a medium voxel grid, and the material of the medium voxel grid is air by default.

[0040] The transport medium can be air, and the target tissue region representing the target object can be specified. For example, the target tissue region can be tumor tissue such as the prostate, liver, lung, or brain. CT images are obtained by scanning the target region of the target object using a tomographic imaging device, and a tissue voxel mesh is constructed based on the CT values ​​of the CT images.

[0041] In a pencil beam ray tracing process, at each determined intersection point, the thread calculates the physical quantities of the path segment leading to that intersection point based on the current AABB material and updates the velocity vector accordingly. The pencil beam then enters the next AABB to iterate this process until energy is exhausted or it exits the computational domain. The pencil beam energy is determined based on the velocity vector. During this process, data is recorded for each pencil beam, such as the indices of tissue voxels traversed within the tissue voxel mesh, the path length from the exit point to each intersection point, and the water-equivalent path length corresponding to the path segment.

[0042] After all threads in the CUDA kernel function for beam tracking have finished executing and exited, the program calls the CUDA kernel function to calculate the dose kernel parameters. In this CUDA kernel function, each GPU thread is assigned to process a pencil beam and is responsible for calculating the dose kernel parameters at all intersection points on its pencil beam trajectory. Specifically, these parameters may include: the standard deviation of the transverse Gaussian broadening caused by multiple Coulomb scattering, the standard deviation of the transverse Gaussian broadening caused by kernel interactions, and the kernel interaction weights. The calculation of these dose kernel parameters is based on data recorded during beam tracking.

[0043] After all threads in the CUDA kernel function responsible for dose kernel parameter calculation have finished executing and exited, the program calls the set of voxels of interest to determine the CUDA kernel function. The set of voxels of interest includes all voxels of interest, which are used for dose calculation. Tissue voxels outside the set of voxels of interest have dose values ​​considered zero and do not need to participate in dose calculation.

[0044] Based on the tissue voxel index obtained from ray tracing, the set of tissue voxels that were traversed was determined. The target distance was formed by calculating the nearest distance from the spatial location of each tissue voxel in the tissue voxel grid to the tissue voxel traversed by the pencil beam. Based on the target distance, tissue voxels whose target distance was less than the field dose cutoff radius were selected and identified as voxels of interest.

[0045] The field dose cutoff radius reflects the upper limit of the action distance at which all pencil beams in a field significantly contribute to the dose distribution in the lateral direction. Specifically, the maximum dose cutoff radius of a single pencil beam is the maximum of its dose cutoff radii at all water equivalent depths. The field dose cutoff radius is determined from the maximum dose cutoff radii of multiple pencil beams. Therefore, the field dose cutoff radius is the maximum of the maximum dose cutoff radii of all pencil beams in the current field.

[0046] A radiation field refers to the set of all pencil beams emitted by a radiotherapy device at a specific gantry angle. A corresponding set of voxels of interest is generated for each radiation field.

[0047] The voxels of interest are defined as tissue voxels within the boundary extension region of all direct pencil beam irradiation or pencil beam lateral scattering ranges under the current field. The boundary extension region is used to fully account for the effects of proton lateral scattering.

[0048] After the CUDA kernel function responsible for determining the voxel of interest has finished executing and all threads have exited, the program immediately calls the CUDA kernel function for the dose calculation task. Within this CUDA kernel function, each GPU thread is assigned a dose calculation subtask to calculate one voxel of interest. The dose calculation subtask completes the calculation by traversing all pencil beams within the current field and summing the voxel deposition dose of each pencil beam for that voxel of interest.

[0049] The calculation process for the voxel deposition dose of a single pencil beam on any voxel of interest is as follows: First, a perpendicular line is drawn from the center of the voxel of interest to the pencil beam trajectory to obtain the foot of the perpendicular; then, based on the position of the foot of the perpendicular, the dose kernel parameters of the pencil beam at multiple intersection points are interpolated to obtain the dose kernel parameters corresponding to the position of the foot of the perpendicular; finally, the dose kernel parameters at the position of the foot of the perpendicular are processed by the dose convolution superposition formula to calculate the voxel deposition dose of the pencil beam on this voxel of interest.

[0050] Based on the doses of multiple pencil beams deposited on a single voxel of interest, the voxel deposition dose corresponding to that voxel is obtained. Based on the individual voxel deposition doses of multiple voxels of interest, the dose distribution results related to a specified tissue are obtained.

[0051] Optionally, since the field dose cutoff radius represents the upper limit of the action distance of all pencil beams in a field that significantly contribute to the dose distribution in the lateral direction, the voxels of interest are selected based on the field dose cutoff radius in the tissue voxel grid. Based on this limitation of the region of interest, the dose superposition calculation is strictly constrained to the area directly irradiated by the pencil beam and its surrounding boundary extension area smaller than the field dose cutoff radius, thereby significantly reducing redundant calculations and significantly improving the dose calculation efficiency. In addition, a voxel-level parallel strategy is adopted for the dose calculation task of multiple voxels of interest. By making full use of the large-scale parallel architecture of the graphics processor, the dose calculation task is decomposed into multiple independently executable dose calculation subtasks, and a large number of threads are scheduled to calculate the voxel deposition dose of the voxels of interest concurrently, further improving the dose calculation efficiency.

[0052] Optionally, based on the field dose cutoff radius, the voxel of interest is determined from the tissue voxels in the tissue voxel grid, including: for any one of the multiple tissue voxels, calculating the distance between the tissue voxel and each of the multiple candidate voxels to obtain multiple distance values, wherein the candidate voxels represent the tissue voxels traversed by the pencil beam; determining the minimum distance value from the multiple distance values; and if the minimum distance value is less than the field dose cutoff radius, determining both the tissue voxel and the candidate voxel corresponding to the minimum distance value as the voxel of interest.

[0053] All tissue voxels traversed by the pen-shaped beam are marked as candidate voxels. The Jump Flooding Algorithm (JFA) is used to calculate the minimum distance value corresponding to each of the multiple tissue voxels in parallel. This process can be as follows: for any tissue voxel and the multiple candidate voxels, multiple distance values ​​are obtained, and the minimum distance value is determined from the multiple distance values.

[0054] When the minimum distance value is less than the cutoff radius of the radiation field dose, this tissue voxel and the candidate voxel corresponding to the minimum distance value are all identified as voxels of interest.

[0055] If the minimum distance value is greater than or equal to the cutoff radius of the radiation field dose, the dose calculation for this tissue voxel will be ignored in subsequent calculations.

[0056] Optionally, a jump flooding algorithm combined with the field dose cutoff radius can be used to screen voxels of interest. The jump flooding algorithm has parallel characteristics and is suitable for GPU implementation. The accuracy of the jump flooding algorithm in calculating distance is sufficient to meet the voxel screening requirements. In addition, the field dose cutoff radius fully considers the influence of proton transverse scattering effect, thereby significantly reducing computational complexity, effectively reducing redundant calculations, and improving the efficiency of voxel screening while ensuring accuracy.

[0057] Optionally, the pencil beam travels in a magnetic field; wherein, calculating the individual pencil beam trajectories of multiple pencil beams includes: using multiple threads to calculate the individual pencil beam trajectories of multiple pencil beams in parallel, wherein any one of the multiple threads is used to perform the following operations: if the velocity value of the pencil beam at the current intersection point is greater than a preset velocity threshold, determine the point magnetic field vector of the current intersection point where the pencil beam is located from a three-dimensional magnetic field distribution map, the three-dimensional magnetic field distribution map representing the spatial magnetic field generated by the magnetic resonance device; determine the upper limit of the deflection angle of the pencil beam at the current intersection point based on the point magnetic field vector, the velocity vector of the pencil beam at the current intersection point, and the straight-line travel distance of the pencil beam within the voxel, the velocity vector including the velocity value; process the point magnetic field vector and the velocity vector according to an intersection point calculation algorithm matching the upper limit of the deflection angle to obtain the next intersection point, the pencil beam trajectory includes multiple intersection points; generate the pencil beam trajectory based on the multiple intersection points.

[0058] The ray tracing task is executed in parallel using multiple threads, with each thread responsible for the ray tracing task of one pen beam, resulting in the pen beam trajectories of multiple pen beams.

[0059] For any thread among multiple threads, first determine whether the pen beam energy is exhausted based on the velocity value of the pen beam at the current intersection point and a preset velocity threshold. If the energy is exhausted, the pen beam trajectory calculation terminates, and the thread exits the calculation. If the velocity value is greater than the preset velocity threshold, it means that the energy is not exhausted. If the energy is not exhausted, continue executing the ray tracing task to determine the next intersection point.

[0060] For example, the preset speed threshold can be 0.

[0061] The 3D magnetic field distribution map is stored in the GPU memory. The data structure of the 3D magnetic field distribution map is a 3D grid, with each grid node storing a magnetic field vector. The data for the 3D magnetic field distribution map comes from pre-measurements or physical simulations of the magnetic field of the magnetic resonance imaging (MRI) device. When the velocity of the pencil beam at the current intersection point exceeds a preset velocity threshold, the point magnetic field vector at the current intersection point is obtained in real-time from the predefined 3D magnetic field distribution map using a trilinear interpolation algorithm, based on the current intersection point of the pencil beam.

[0062] Figure 3 A schematic diagram illustrating a two-dimensional magnetic field distribution according to an embodiment of the present invention is shown.

[0063] To clearly demonstrate the magnetic field environment used in the embodiments of the present invention, Figure 3 This invention provides a magnetic field distribution based on three-dimensional magnetic field data in the z=0cm plane. The horizontal axis represents the positional information along the x-axis, and the vertical axis represents the positional information along the y-axis. The pixel values ​​of the two-dimensional magnetic field distribution map represent the magnetic field strength. This illustration is for illustrative purposes only and is intended to aid in understanding the macroscopic distribution of the spatial magnetic field; it is not intended for actual calculations. The three-dimensional magnetic field distribution map used in this embodiment of the invention is magnetic field data covering the entire three-dimensional space obtained by physically modeling a 1.5T split-type superconducting MRI coil for clinical applications using modeling software, and is used for actual calculations.

[0064] The point magnetic field vector includes the magnetic field strength and direction at the current intersection point.

[0065] For the current intersection point, a voxel is defined that contains the current intersection point, and the pencil beam enters this voxel along the current velocity direction. The velocity vector at the current intersection point includes the velocity direction. The distance traveled by the pencil beam within the voxel represents the geometric distance the pencil beam traverses within that voxel. Assuming that the proton propagates in a straight line, the straight-line travel distance of the proton within the voxel is calculated using a straight-line intersection algorithm.

[0066] Calculate the maximum radius of curvature of the trajectory of the pencil beam under the point magnetic field vector, based on the point magnetic field vector, the velocity vector of the pencil beam at the current intersection point, and the distance traveled by the pencil beam within the voxel.

[0067] In one embodiment, the upper limit of the deflection angle of the pencil beam at the current intersection point. As shown in formula (1):

[0068] (1).

[0069] in, Characterizes the linear travel distance of the pencil beam within a voxel. The maximum radius of curvature of the trajectory of the pen-shaped beam is characterized by the velocity vector and the point magnetic field vector at the current intersection point. The coefficient 1.1 is a conservative coefficient used to fully account for the errors that may be introduced by the straight-line approximation and to ensure that a certain degree of safety is maintained in the deflection estimation.

[0070] The curvature of the pencil beam trajectory near the current intersection point is evaluated based on the upper limit of the deflection angle. A matching intersection point calculation algorithm is then selected to process the point magnetic field vector and velocity vector of the pencil beam, obtaining the next intersection point. The velocity vector of the proton at the next intersection point is updated, and multiple intersection points are obtained sequentially. The pencil beam trajectory calculation terminates when the calculated intersection point is not in the magnetic field region, or when the pencil beam energy is exhausted, and the thread exits the calculation.

[0071] Pencil beam trajectories are generated based on multiple intersection points located within the magnetic field region.

[0072] Optionally, during ray tracing, a pencil beam parallelization strategy is adopted, with each pencil beam processed by an independent thread. This thread iteratively calculates multiple intersection points based on the pencil beam's velocity vector, point magnetic field vector, and travel distance within the voxel, scheduling a large number of threads to execute concurrently, fully utilizing the GPU's parallel computing capabilities. Furthermore, a dynamic memory resource management strategy is implemented for the 3D magnetic field distribution map. By strictly limiting calls to the 3D magnetic field distribution map to the execution scope of the ray tracing CUDA kernel function and immediately releasing its memory after the CUDA kernel function exits, unnecessary long-term memory occupation is effectively avoided, thereby minimizing the pressure on GPU memory resources from large-scale data and ensuring the smoothness of the overall computation process.

[0073] Optionally, the intersection calculation algorithm is determined based on the following operations: if the upper limit of the deflection angle is less than the preset deflection angle threshold, the intersection calculation algorithm is determined to be a straight line intersection algorithm; if the upper limit of the deflection angle is greater than or equal to the preset deflection angle threshold, the intersection calculation algorithm is determined to be a spiral line intersection algorithm.

[0074] For example, the preset deflection angle threshold can be 0.01 rad (radians).

[0075] If the upper limit of the deflection angle is less than the preset deflection angle threshold, it means that the deflection of the pencil beam within the currently moving voxel is negligible. The trajectory of the pencil beam near the current intersection point is approximately linear. Therefore, the intersection point calculation algorithm is a straight-line intersection algorithm, thus avoiding the costly iterative numerical root-finding process. The pencil beam travels along a straight line within the voxel, and its travel distance is denoted as the straight-line travel distance.

[0076] If the upper limit of the deflection angle is greater than or equal to the preset deflection angle threshold, it means that the deflection of the pen beam within the currently moving voxel is not negligible, and the trajectory of the pen beam near the current intersection point is approximately a spiral trajectory. Therefore, the intersection point calculation algorithm is a spiral intersection algorithm. The pen beam travels along an arc within the voxel, and its travel distance is the arc travel distance.

[0077] For example, helix intersection algorithms can include numerical root-finding algorithms.

[0078] Optionally, the main magnetic field of MRI is a non-uniform magnetic field, which can cause complex deflection of the pencil beam, significantly increasing the difficulty of dose calculation. Therefore, this invention designs an adaptive mechanism based on the assessment of deflection degree in the ray tracing process. This mechanism adaptively switches the corresponding intersection calculation method according to the curvature of the local trajectory of the pencil beam. This design is applicable to both parallel and perpendicular magnetic field configurations, strictly ensuring tracking accuracy in high-curvature critical regions (such as near the Bragg peak), while avoiding unnecessary numerical root calculations in most low-curvature regions, thereby achieving a significant improvement in computational efficiency on a global scale.

[0079] Optionally, the velocity vector of the pencil beam at the current intersection point is determined based on the following operation: the velocity vector of the pencil beam at the current intersection point is determined according to the deflection angle, point magnetic field vector, and velocity vector of the pencil beam at the previous intersection point.

[0080] Update the velocity vector of the pencil beam at the current intersection point based on the deflection angle, point magnetic field vector, and velocity vector of the pencil beam at the previous intersection point.

[0081] The deflection angle at the previous intersection point is related to the intersection calculation algorithm used to calculate the current intersection point. Specifically, the upper limit of the deflection angle at the previous intersection point and the preset deflection angle threshold are compared and a matching intersection calculation algorithm is selected. The deflection angle of the pencil beam at the previous intersection point is calculated using the travel distance matched with the intersection calculation algorithm. This intersection calculation algorithm is used to process the point magnetic field vector and velocity vector of the previous intersection point to obtain the current intersection point.

[0082] Depending on the intersection calculation algorithm, the travel distance can be either a straight line travel distance or an arc travel distance. For example, when the intersection calculation algorithm is a straight line intersection algorithm, the travel distance is the straight line travel distance of the proton within the voxel; when the intersection calculation algorithm is a spiral intersection algorithm, the travel distance is the arc travel distance of the proton within the voxel.

[0083] In one embodiment, the deflection angle of the pencil beam at the previous intersection point As shown in formula (2):

[0084] (2).

[0085] in, Characterizes the distance traveled by the pen-shaped beam within the voxel associated with the previous intersection point.

[0086] During the iteration process, an intersection point is determined each time, and the velocity vector of the pen beam at this intersection point needs to be updated.

[0087] Optionally, based on the dose kernel parameters of each of the multiple intersection points, multiple threads are invoked to execute dose calculation tasks for multiple voxels of interest in parallel to obtain the voxel deposition dose of the voxels of interest. This includes: using any one of the multiple threads to perform the following operations for any pencil beam: determining the radial distance from the voxel center point to the foot of the perpendicular from the voxel center point to the pencil beam trajectory, and the path length of the pencil beam from the exit point to the foot of the perpendicular; when the path length is less than the range of the pencil beam and the radial distance is less than the dose cutoff radius of the pencil beam, using an interpolation algorithm to process the path length and the dose kernel parameters of the pencil beam at the multiple intersection points to obtain the dose kernel parameters corresponding to the foot of the perpendicular; obtaining the dose of the pencil beam deposited on the voxel of interest based on the dose of each of the multiple pencil beams deposited on the voxel of interest; and obtaining the voxel deposition dose of each of the multiple voxels of interest in parallel based on multiple threads.

[0088] For any one of the multiple threads, the task of calculating the dose of the voxel of interest by this thread includes: traversing all pencil beams in the current field and obtaining the voxel deposition dose corresponding to the voxel of interest based on the sum of the doses deposited by each pencil beam in this voxel of interest.

[0089] For example, for any pen beam P among multiple pen beams, draw a perpendicular line from the voxel center point of the voxel of interest to the trajectory of pen beam P to obtain the foot of the perpendicular to the trajectory of pen beam P, and then calculate the distance between the voxel center point and the foot of the perpendicular, which is the radial distance.

[0090] The path length is determined based on the length of the pen-shaped beam P from the exit point to the foot of the perpendicular.

[0091] The preset range threshold represents the upper limit of the pen beam's range before its energy is depleted.

[0092] The dose cutoff radius of a pencil beam characterizes the upper limit of the distance at which the pencil beam makes a significant contribution to the dose distribution in the lateral direction.

[0093] When the path length is less than the preset range threshold of the pencil beam P and the radial distance is less than the dose cutoff radius of the pencil beam P, the path length and the dose kernel parameters of the pencil beam at multiple intersection points are processed by the interpolation algorithm to obtain the dose of the pencil beam P at the voxel of interest.

[0094] Specifically, the dose kernel parameter of the foot of the voxel of interest depends on the path length from the exit point of the pencil beam along its trajectory to the foot. Based on this, for the foot of the voxel center on the pencil beam trajectory, the dose kernel parameter corresponding to the foot is obtained by interpolating the dose kernel parameters of its multiple neighboring intersection points. Then, the dose kernel parameter of the foot is processed using the dose convolution superposition formula to obtain the dose of the pencil beam P at the voxel of interest.

[0095] The dose of each pencil beam deposited on the voxel of interest under the current field is calculated sequentially, and the dose contributions of all pencil beams are summed to obtain the voxel deposition dose of the voxel of interest.

[0096] Based on multiple threads, the voxel deposition dose of each voxel of interest is obtained in parallel.

[0097] Optionally, a voxel parallel strategy is adopted in the dose calculation task, with each GPU thread responsible for calculating the voxel deposition dose of a voxel of interest, thereby achieving large-scale parallel acceleration. This invention integrates two key optimizations in the dose calculation task to improve efficiency. First, pre-screening is implemented: for each pencil beam, based on the radial distance from the voxel center to the foot of the pencil beam trajectory and the path length from the exit point to the foot, it is determined whether its dose contribution is negligible; if so, subsequent calculations are immediately skipped to avoid redundant computation. Second, data reuse is achieved through dose kernel parameter interpolation: for pencil beams with dose contributions, based on the path length from the exit point to the foot, a linear interpolation algorithm is used to directly obtain the dose kernel parameters at the foot from the dose kernel parameters of multiple intersection points, avoiding the expensive recalculation of dose kernel parameters at the foot position, further improving computational efficiency.

[0098] Optionally, the dose kernel parameters include the Gaussian standard deviation of nuclear interactions, the Gaussian standard deviation of multiple Coulomb scattering, and the nuclear interaction weights. The dose kernel parameters for each pencil beam at multiple intersection points are obtained based on the following operation: for any intersection point, the Gaussian standard deviation of multiple Coulomb scattering is determined according to the target beam size and the first beam size. The target beam size characterizes the lateral broadening of the pencil beam as it travels through voxels, and the first beam size characterizes the lateral broadening of the pencil beam in tissue voxels due to multiple Coulomb scattering. The first beam size is a base... The path length of the pencil beam from the initial tissue intersection point to the current intersection point is determined, where the initial tissue intersection point represents the first intersection point between the pencil beam and the tissue voxel mesh. The Gaussian standard deviation of the nuclear interaction is determined based on the Gaussian standard deviation of multiple Coulomb scattering and the second beam spot size, where the second beam spot size represents the lateral broadening of the pencil beam in the tissue voxel due to nuclear interaction, and the second beam spot size is determined based on the path length of the pencil beam from the initial tissue intersection point to the current intersection point. The nuclear interaction weight is determined based on the path length of the pencil beam from the initial tissue intersection point to the current intersection point.

[0099] The target beam spot size characterizes the lateral widening caused by the pen-shaped beam traveling in tissue voxels or medium voxels.

[0100] The path length of the pencil beam from the initial tissue intersection point to the current intersection point is the geometric distance from the initial tissue intersection point to the current intersection point. The square of the Gaussian standard deviation of the multiple Coulomb scattering interactions is determined by the sum of the squares of the target beam size and the squares of the first beam size.

[0101] The square of the Gaussian standard deviation of the nuclear interaction is determined by summing the square of the Gaussian standard deviation of the multiple Coulomb scattering interactions and the square of the second beam spot size.

[0102] Based on the path length of the foot of the perpendicular corresponding to the voxel of interest, the Gaussian standard deviation of the nuclear interaction, the Gaussian standard deviation of the multiple Coulomb scattering, and the nuclear interaction weights of each of the multiple intersection points are interpolated to obtain the Gaussian standard deviation of the nuclear interaction, the Gaussian standard deviation of the multiple Coulomb scattering, and the nuclear interaction weights of the foot of the perpendicular. Then, the Gaussian distribution component of the nuclear interaction is determined based on the Gaussian standard deviation of the nuclear interaction, and the Gaussian distribution component of the multiple Coulomb scattering is determined based on the Gaussian standard deviation of the multiple Coulomb scattering.

[0103] Based on the Gaussian distribution components of nuclear interactions, the Gaussian distribution components of multiple Coulomb scattering, and the nuclear interaction weights, the perpendicular double Gaussian dose kernel function is obtained.

[0104] The nuclear interaction weight is the weight of the total reaction in which nuclear interactions occur.

[0105] In one embodiment, the bigaussian dose kernel function at the foot of the perpendicular As shown in formula (3):

[0106] (3).

[0107] in, Characterizing nuclear interaction weights, The Gaussian distribution component characterizing multiple Coulomb scattering effects. The Gaussian distribution component characterizing nuclear interactions, r characterizing the radial distance from the center of the voxel of interest to the foot of the pendulum trajectory. It represents the path length from the exit point of the pencil beam to the foot of the perpendicular.

[0108] In one embodiment, the pencil beam is used to dose the voxel of interest. As shown in formula (4):

[0109] (4).

[0110] in, Characterizing the integral depth dose, This is a double Gaussian dose kernel function.

[0111] Optionally, the target beam size is obtained based on the following operations: the initial optical parameters are looked up from the optical parameter lookup table based on the initial energy of the pencil beam, the strength of the main magnetic field, and the direction of the magnetic field configuration; the target beam size is determined based on the initial optical parameters and the path length of the pencil beam from the exit point to the current intersection point.

[0112] Under the influence of a non-uniform edge magnetic field in MRI, the proton beam will exhibit focusing: when the magnetic field is parallel to the incident direction of the radiation field, the beam mainly rotates around the axis and is accompanied by a slight focusing effect; when the magnetic field is perpendicular to the incident direction of the radiation field, the pencil beam will not only deflect laterally, but will also be compressed in the direction parallel to the magnetic field.

[0113] In MRPT, the target object is placed at the center of the MRI main magnetic field (imaging field). The magnetic field in this region has high homogeneity, so the magnetic focusing effect can be ignored when the pencil beam propagates within the target object. The magnetic focusing effect mainly occurs when the pencil beam passes through a non-uniform edge magnetic field region. This process changes the optical parameters of the beam, such as beam spot size, angular divergence, and spatial-angular covariance.

[0114] To reproduce the specific physical scenario of a pencil beam simultaneously traversing both the imaging field and the edge field in real clinical settings, and to verify the applicability of the algorithm under these complex conditions, this invention positions the pencil beam phase space source 1.5 m from the isocenter point. This setup effectively triggers the deflection and focusing effect of the proton beam in the magnetic field, thereby fully verifying the accuracy and applicability of the dose calculation method in a real treatment environment.

[0115] The objective of this invention is to accurately calculate the dose distribution of a pencil beam in a specified tissue of a target object. Therefore, it is not necessary to model the intermediate process of magnetic focusing; only the target beam spot size of the pencil beam needs to be obtained to calculate the dose kernel parameters.

[0116] Each magnetic field has a fixed main magnetic field strength and direction (magnetic field configuration direction). The initial optical parameters are obtained by consulting an initial optical parameter lookup table based on the main magnetic field strength, the initial energy of the pencil beam, and the magnetic field configuration direction. This lookup table is a pre-generated Monte Carlo simulation lookup table. It takes the magnetic field configuration direction (parallel / perpendicular), the main magnetic field strength, and the initial energy of the pencil beam as input, and outputs the corresponding initial optical parameters for the pencil beam. The initial energy of the pencil beam is the energy emitted from the beam.

[0117] The initial optical parameters include the initial beam size, angular divergence, and spatial-angular covariance. Then, based on the initial beam size, angular divergence, spatial-angular covariance, and path length of the intersection point, the target beam size at the current position is calculated.

[0118] The initial beam spot size is the initial spatial variance, which reflects the lateral size of the beam at the incident surface.

[0119] Angular divergence is the initial angular variance, reflecting the degree of beam divergence.

[0120] Spatial-angular covariance is used to describe the relationship between the lateral position of the beam and the exit direction.

[0121] In one embodiment, the target beam spot size As shown in formula (5):

[0122] (5).

[0123] in, Characterized by the initial beam spot size, For space-angle covariance, Let be the path length of the pen-shaped beam from the exit point to the current intersection point. The representation perspective is divergent.

[0124] The path length from the exit point to the current intersection point is the sum of the path length from the exit point to the initial tissue intersection point and the path length from the initial tissue intersection point to the current intersection point.

[0125] Optionally, this invention directly performs high-precision dose calculation based on the real and complete three-dimensional magnetic field distribution of clinical split MRI, and simultaneously supports parallel and perpendicular magnetic field configurations. This overcomes the dependence of existing dose algorithms on an ideal uniform perpendicular magnetic field, thus possessing broader clinical applicability. Furthermore, this invention accurately characterizes the magnetic focusing effect caused by the edge field by calculating the target beam spot size using an initial optical parameter lookup table. The focusing effect of the edge magnetic field is here equivalent to the modulation of the initial optical parameters, thereby effectively characterizing the influence of the edge magnetic field on pencil beam transmission while maintaining algorithm simplicity.

[0126] Optionally, the field dose cutoff radius is obtained based on the following operations: For any pencil beam among multiple pencil beams in a single field, the pencil beam trajectory is sampled according to a preset sampling rule to obtain multiple sampling points; based on the lateral dose distribution of the sampling points, the dose cutoff radius corresponding to the sampling points is obtained using a lateral dose integral function, where the lateral dose distribution is obtained by interpolating the dose kernel parameters of multiple intersection points on the pencil beam trajectory; the maximum dose cutoff radius among the dose cutoff radii of the multiple sampling points is determined as the dose cutoff radius of the pencil beam; and the maximum dose cutoff radius among the dose cutoff radii of the multiple pencil beams is determined as the field dose cutoff radius.

[0127] The preset sampling rule is to select 5 equidistant points for sampling within 30 mm before the end of the pencil beam trajectory (where the proton lateral scattering broadening is most significant).

[0128] For each sampling point, the dose kernel parameters of multiple intersection points near the sampling point are processed by a linear interpolation algorithm to obtain the dose kernel parameters of the sampling point, thereby obtaining the lateral dose distribution.

[0129] In one embodiment, the transverse dose integral function is shown in equation (6):

[0130] (6);

[0131] in, The path length from the pencil beam exit point to the current sampling point is used as the characteristic parameter corresponding to the sampling point. In formula (6), each sampling point has its corresponding fixed value. Characterizes the radial distance relative to the pencil beam trajectory; Characterized in path length Lateral dose distribution at point 2; The minimum radial distance that satisfies formula (6), i.e., the path length, is represented by this formula. The dose cutoff radius at the sampling point; Characterizing the radial differential; The cutoff threshold is used to represent the proportion of the dose integral that needs to be retained; for example, it can be set... =0.99. Formula (6) means: select a minimum radial distance. This ensures that the proportion of the cumulative deposition dose within this radial distance to the total lateral deposition dose is not less than the cutoff threshold. .

[0132] For example, when the truncation level threshold is reached When the path length is At the sampling point, the lateral cumulative dose within the dose cutoff radius is not less than 99% of the total lateral dose integral at that sampling point. The remaining approximately 1% lateral dose contribution can be ignored.

[0133] The maximum dose cutoff radius among the dose cutoff radii of multiple sampling points is determined as the dose cutoff radius of the pencil beam, and the maximum dose cutoff radius among the dose cutoff radii of multiple pencil beams is determined as the field dose cutoff radius.

[0134] Figure 4 A flowchart of a proton dose distribution determination method based on a graphics processor according to another embodiment of the present invention is shown.

[0135] like Figure 4 As shown, the process includes steps S410 to S480.

[0136] In step S410, the GPU loads the target object image data, the proton beam phase space source, the optical parameter lookup table, and the three-dimensional magnetic field distribution map, allocates memory, and transfers the required data from the host memory to the device memory.

[0137] In step S420, pen-beam ray tracing is performed, with each thread responsible for calculating the trajectory of one pen-beam. Specifically, step S420 includes steps S421 to S428.

[0138] In step S421, it is determined whether the speed of the pen beam is greater than zero. If so, this thread exits the calculation; otherwise, proceed to step S422.

[0139] In step S422, the point magnetic field vector at the current intersection point of the pencil beam is obtained from the three-dimensional magnetic field distribution map.

[0140] In step S423, the upper limit of the deflection angle of the pencil beam at the current intersection point is determined based on the point magnetic field vector, the velocity vector, and the linear travel distance of the pencil beam within the voxel.

[0141] In step S424, it is determined whether the upper limit of the deflection angle is less than 0.01 radians. If so, proceed to step S425; otherwise, proceed to step S426.

[0142] In step S425, the point magnetic field vector and velocity vector are processed using a straight line intersection algorithm to obtain the next intersection point, and the index of the voxel that is crossed is recorded.

[0143] In step S426, the point magnetic field vector and velocity vector are processed using the spiral intersection algorithm to obtain the next intersection point, and the index of the voxel that is crossed is recorded.

[0144] In step S427, determine whether the next phase intersection point is still within the magnetic field region; if so, proceed to step S428; otherwise, this thread exits the calculation.

[0145] In step S428, the velocity vector of the proton at the next phase intersection point is updated, and the process proceeds to step S421.

[0146] After all threads simulating pencil beam trajectories exit the calculation, proceed to step S430.

[0147] In step S430, the dose kernel parameters of each of the multiple intersection points are obtained.

[0148] In step S440, the dose cutoff radius of each pencil beam and the dose cutoff radius of the radiation field are determined.

[0149] In step S450, the voxel of interest is determined.

[0150] In step S460, a dose calculation task is performed, with each thread responsible for calculating the deposition dose of the voxel of interest. Specifically, step S460 includes steps S461 to S467.

[0151] In step S461, the thread obtains the voxel center point of the voxel it is responsible for and initializes the pen beam index value i = 1.

[0152] In step S462, determine whether all pen-shaped beams in the current field of view have been calculated: if so, the thread exits the calculation; otherwise, proceed to step S463.

[0153] In step S463, the radial distance from the voxel center point to the foot of the perpendicular and the path length of the pencil beam from the exit point to the foot of the perpendicular are calculated.

[0154] In step S464, it is determined whether the path length is less than the preset range threshold of the pencil beam and whether the radial distance is less than the dose cutoff radius of the pencil beam. If both are satisfied, proceed to step S465; otherwise, increment the pencil beam index i by 1 and proceed to step S462.

[0155] In step S465, based on the path length, an interpolation algorithm is used to process the path length and the dose kernel parameters of the pencil beam at multiple intersection points to obtain the dose of the pencil beam deposited on the voxel of interest.

[0156] In step S466, the dose contribution of the i-th pencil beam to the voxel of interest is calculated, and the dose is added to the total dose of the voxel of interest.

[0157] In step S467, the pen-shaped bundle index i is incremented by 1, and the process proceeds to S462.

[0158] After all threads performing dose calculation tasks exit the calculation, proceed to step S470.

[0159] In step S470, determine whether all shooting fields have been processed: if yes, proceed to step S480; otherwise, proceed to step S450.

[0160] In step S480, the dose distribution results are output.

[0161] Figure 5 A schematic diagram showing the dose distribution results of a simulated proton incident on a water model according to an embodiment of the present invention is illustrated.

[0162] like Figure 5 As shown, under a 1.5 T split-type MRI three-dimensional magnetic field, the dose distribution results obtained by proton beam incident on a water phantom using parallel and perpendicular magnetic field configurations, respectively. Figure 5 a and Figure 5 d represents the dose distribution result calculated by the proton dose distribution determination method based on the graphics processor of this invention; Figure 5 b and Figure 5 e represents the dose distribution result based on the Monte Carlo method, serving as a baseline reference; Figure 5 c and Figure 5 f is a pass rate distribution graph showing the dose distribution obtained by the method of this invention and the baseline reference dose distribution after Gamma analysis based on the 2 mm / 2% standard. Specifically, Figure 5 a, Figure 5 b、 Figure 5 c corresponds to a parallel magnetic field configuration. Figure 5 d、 Figure 5 e Figure 5 f corresponds to the vertical magnetic field configuration. Analysis results show that the dose distribution calculated by the method of this invention is highly consistent with the dose distribution results of the Monte Carlo method, verifying that the method of this invention still has high computational accuracy under complex magnetic field environments. Gamma analysis shows that, under a 10% threshold condition, the 2 mm / 2% Gamma pass rate reaches 100% under both magnetic field configurations. In terms of computational efficiency, the method of this invention exhibits a significant advantage: completing... Figure 5 a and Figure 5 The three-dimensional dose calculation shown in d only requires 2887 ms (parallel field) and 3789 ms (vertical field), respectively, which fully meets the efficiency requirements of clinical online adaptive applications.

[0163] Figure 6 A schematic diagram showing the dose distribution results of simulated proton incident on a liver tumor according to an embodiment of the present invention is illustrated.

[0164] like Figure 6 The image shows the dose distribution results obtained when a proton beam is incident on the liver region under a 1.5 T split-type MRI three-dimensional magnetic field, using parallel and perpendicular magnetic field configurations, respectively. Among them, Figure 6 a and Figure 6 d represents the dose distribution result calculated by the proton dose distribution determination method based on the graphics processor of this invention; Figure 6 b and Figure 6 e represents the dose distribution result simulated using the Monte Carlo method, serving as a baseline reference; Figure 6 c and Figure 6 f is a pass rate distribution graph showing the dose distribution obtained by the method of this invention and the baseline reference dose distribution after Gamma analysis based on the 2 mm / 2% standard. Specifically, Figure 6 a, Figure 6 b、 Figure 6 c corresponds to a parallel magnetic field configuration. Figure 6 d、 Figure 6 e Figure 6 f corresponds to the vertical magnetic field configuration. In order to ensure that the proton dose can cover the tumor target area in a magnetic field environment, the original treatment plan (i.e., the plan formulated under non-magnetic field conditions) was adjusted.

[0165] To compensate for the influence of the magnetic field on the pencil beam trajectory and ensure that protons fall within the target area, the treatment plan developed under the original no-magnetic-field conditions was adjusted. In the parallel magnetic field configuration, the point distribution on the virtual source plane was rotated by 18° to correct the beam rotation under the influence of the edge magnetic field. In the vertical magnetic field configuration, the gantry of the two radiation fields was rotated by 25° and 26° respectively based on the original gantry angle, and the isocenter positions were translated by 80 mm and 35 mm respectively in the longitudinal direction to compensate for the beam deflection caused by the magnetic field.

[0166] The analysis results show that the dose distribution calculated by the method of this invention is highly consistent with the dose distribution simulated by the Monte Carlo method, verifying that the method of this invention still has high computational accuracy under complex magnetic field environments. Gamma analysis shows that, under the 10% dose threshold condition, the throughput in parallel field and vertical field configurations reaches 99.88% and 99.99%, respectively. In terms of computational efficiency, the method of this invention exhibits significant advantages: the three-dimensional dose calculation for this liver tumor case only requires 3439 ms (parallel field) and 3658 ms (vertical field), fully meeting the efficiency requirements of clinical online adaptive applications.

[0167] Based on the above-described proton dose distribution determination method based on a graphics processor, this invention also provides a proton dose distribution determination device based on a graphics processor. The following will be combined with... Figure 7 The device is described in detail.

[0168] Figure 7 A structural block diagram of a proton dose distribution determination device based on a graphics processor according to an embodiment of the present invention is shown.

[0169] like Figure 7 As shown, the proton dose distribution determination device 700 based on a graphics processor in this embodiment includes a trajectory calculation module 710, an acquisition module 720, a first determination module 730, a calling module 740, and a second determination module 750.

[0170] The trajectory calculation module 710 is used to calculate the trajectory of each of the multiple pencil beams, wherein the intersection point in the pencil beam trajectory represents the intersection of the pencil beam with the medium voxel grid representing the transmission medium, or the tissue voxel grid representing a specified tissue.

[0171] The acquisition module 720 is used to acquire the dose kernel parameters of the pencil beam at multiple intersection points.

[0172] The first determining module 730 is used to determine the voxel of interest from tissue voxels in the tissue voxel grid according to the field dose cutoff radius, wherein the field dose cutoff radius characterizes the upper limit of the action distance of all pencil beams under a field that make a significant contribution to the dose distribution in the lateral direction.

[0173] Module 740 is invoked to call multiple threads to perform dose calculation tasks on multiple voxels of interest in parallel based on the dose kernel parameters of each of the multiple intersection points, so as to obtain the voxel deposition dose of the voxels of interest.

[0174] The second determination module 750 is used to determine the dose distribution results related to a specified tissue based on the voxel deposition doses of the multiple voxels of interest.

[0175] Optionally, the first determining module 730 includes a first determining submodule, a second determining submodule, and a third determining submodule.

[0176] The first determination submodule is used to calculate the distance between any one of the multiple tissue voxels and each of the multiple candidate voxels, and obtain multiple distance values, wherein the candidate voxels represent tissue voxels that are traversed by the pen-shaped beam.

[0177] The second determination submodule is used to determine the minimum distance value from multiple distance values.

[0178] The third determination submodule is used to determine both the tissue voxel and the candidate voxel corresponding to the minimum distance value as the voxel of interest when the minimum distance value is less than the cutoff radius of the radiation field dose.

[0179] Optionally, any multiple modules among the trajectory calculation module 710, acquisition module 720, first determination module 730, invocation module 740, and second determination module 750 can be combined into one module, or any one of these modules can be split into multiple modules. Alternatively, at least part of the functionality of one or more of these modules can be combined with at least part of the functionality of other modules and implemented in one module. Optionally, at least one of the trajectory calculation module 710, acquisition module 720, first determination module 730, invocation module 740, and second determination module 750 can be at least partially implemented as hardware circuitry, such as a field-programmable gate array (FPGA), programmable logic array (PLA), system-on-a-chip, system-on-a-substrate, system-on-package, application-specific integrated circuit (ASIC), or any other reasonable means of integrating or packaging circuitry, or implemented in software, hardware, or firmware, or in any appropriate combination of any of these three implementation methods. Alternatively, at least one of the trajectory calculation module 710, acquisition module 720, first determination module 730, calling module 740 and second determination module 750 can be at least partially implemented as a computer program module, which can perform corresponding functions when the computer program module is run.

[0180] Figure 8 A block diagram of an electronic device suitable for implementing a graphics processor-based proton dose distribution determination method according to an embodiment of the present invention is shown.

[0181] Figure 8 The electronic device shown is merely an example and should not be construed as limiting the functionality and scope of use of the embodiments of the present invention.

[0182] like Figure 8 As shown, a computer electronic device 800 according to an embodiment of the present invention includes a processor 801, which can perform various appropriate actions and processes according to a program stored in a ROM 802 (read-only memory) or a program loaded from a storage portion 808 into a RAM 803 (random access memory). The processor 801 may include, for example, a general-purpose microprocessor (e.g., a CPU), an instruction set processor and / or an associated chipset and / or a special-purpose microprocessor (e.g., an application-specific integrated circuit (ASIC)), etc. The processor 801 may also include onboard memory for caching purposes. The processor 801 may include a single processing unit or multiple processing units for performing different actions of the method flow according to an embodiment of the present invention.

[0183] RAM 803 stores various programs and data required for the operation of electronic device 800. Processor 801, ROM 802, and RAM 803 are interconnected via bus 804. Processor 801 executes various operations of the method flow according to embodiments of the present invention by executing programs in ROM 802 and / or RAM 803. It should be noted that programs may also be stored in one or more memories other than ROM 802 and RAM 803. Processor 801 may also execute various operations of the method flow according to embodiments of the present invention by executing programs stored in one or more memories.

[0184] Optionally, the electronic device 800 may also include an input / output (I / O) interface 805, which is also connected to the bus 804. The electronic device 800 may also include one or more of the following components connected to the input / output (I / O) interface 805: an input section 806 including a keyboard, mouse, etc.; an output section 807 including a cathode ray tube (CRT), liquid crystal display (LCD), etc., and speakers, etc.; a storage section 808 including a hard disk, etc.; and a communication section 809 including a network interface card such as a LAN card, modem, etc. The communication section 809 performs communication processing via a network such as the Internet. A drive 810 is also connected to the input / output (I / O) interface 805 as needed. A removable medium 811, such as a disk, optical disk, magneto-optical disk, semiconductor memory, etc., is installed on the drive 810 as needed so that computer programs read from it can be installed into the storage section 808 as needed.

[0185] Optionally, the method flow according to embodiments of the present invention can be implemented as a computer software program. For example, embodiments of the present invention include a computer program product comprising a computer program carried on a computer-readable storage medium, the computer program containing program code for performing the method shown in the flowchart. In such embodiments, the computer program can be downloaded and installed from a network via communication section 809, and / or installed from removable medium 811. When the computer program is executed by processor 801, it performs the functions defined in the system of embodiments of the present invention. Optionally, the systems, devices, apparatuses, modules, units, etc., described above can be implemented by computer program modules.

[0186] The present invention also provides a computer-readable storage medium, which may be included in the device / apparatus / system described in the above embodiments; or it may exist independently and not assembled into the device / apparatus / system. The computer-readable storage medium carries one or more programs, which, when executed, implement the graphics processor-based proton dose distribution determination method according to embodiments of the present invention.

[0187] Optionally, the computer-readable storage medium can be a non-volatile computer-readable storage medium. Examples include, but are not limited to: portable computer disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), portable compact disk read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination thereof. In this invention, the computer-readable storage medium can be any tangible medium containing or storing a program that can be used by or in conjunction with an instruction execution system, apparatus, or device.

[0188] For example, optionally, the computer-readable storage medium may include the ROM 802 and / or RAM 803 described above and / or one or more memories other than ROM 802 and RAM 803.

[0189] Embodiments of the present invention also include a computer program product comprising a computer program containing program code for performing the methods provided in the embodiments of the present invention. When the computer program product is run on an electronic device, the program code is used to enable the electronic device to implement the graphics processor-based proton dose distribution determination method provided in the embodiments of the present invention.

[0190] When the computer program is executed by the processor 801, it performs the functions defined in the system / apparatus of this embodiment of the invention. Optionally, the systems, apparatuses, modules, units, etc., described above can be implemented by computer program modules.

[0191] In one embodiment, the computer program may rely on a tangible storage medium such as an optical storage device or a magnetic storage device. In another embodiment, the computer program may also be transmitted and distributed in the form of signals over a network medium, and may be downloaded and installed via the communication section 809, and / or installed from a removable medium 811. The program code contained in the computer program can be transmitted using any suitable network medium, including but not limited to: wireless, wired, etc., or any suitable combination thereof.

[0192] Optionally, program code for executing the computer programs provided in the embodiments of the present invention can be written in any combination of one or more programming languages. Specifically, these computational programs can be implemented using high-level procedural and / or object-oriented programming languages, and / or assembly / machine languages. Programming languages ​​include, but are not limited to, languages ​​such as Java, C++, Python, "C", or similar programming languages. The program code can be executed entirely on the user's computing device, partially on the user's device, partially on a remote computing device, or entirely on a remote computing device or server. In cases involving remote computing devices, the remote computing device can be connected to the user's computing device via any type of network, including a local area network (LAN) or a wide area network (WAN), or it can be connected to an external computing device (e.g., via the Internet using an Internet service provider).

[0193] The flowcharts and block diagrams in the accompanying drawings illustrate the architecture, functionality, and operation of possible implementations of systems, methods, and computer program products according to various embodiments of the present invention. In this regard, each block in a flowchart or block diagram may represent a module, segment, or portion of code containing one or more executable instructions for implementing a specified logical function. It should also be noted that in some alternative implementations, the functions indicated in the blocks may occur in a different order than those indicated in the drawings. For example, two consecutively indicated blocks may actually be executed substantially in parallel, and they may sometimes be executed in reverse order, depending on the functions involved. It should also be noted that each block in a block diagram or flowchart, and combinations of blocks in a block diagram or flowchart, may be implemented using a dedicated hardware-based system that performs the specified function or operation, or using a combination of dedicated hardware and computer instructions. Those skilled in the art will understand that the features described in the various embodiments of the present invention can be combined and / or combined in various ways, even if such combinations or combinations are not explicitly described in the present invention. In particular, the features described in the various embodiments of the present invention can be combined and / or combined in various ways without departing from the spirit and teachings of the present invention. All such combinations and / or pairings fall within the scope of this invention.

[0194] The embodiments of the present invention have been described above. However, these embodiments are merely illustrative and not intended to limit the scope of the invention. Although various embodiments have been described above, this does not mean that the measures in the various embodiments cannot be used advantageously in combination. Various substitutions and modifications can be made by those skilled in the art without departing from the scope of the invention, and all such substitutions and modifications should fall within the scope of the invention.

Claims

1. A method for determining proton dose distribution based on a graphics processor, characterized in that, include: Calculate the pencil beam trajectory of each of the multiple pencil beams, wherein the intersection points in the pencil beam trajectory include the intersection points where the pencil beam intersects with the medium voxel grid of the transmission medium, and the intersection points where the pencil beam intersects with the tissue voxel grid of a specified tissue, wherein the transmission medium is air; Obtain the dose kernel parameters of the pencil beam at multiple intersection points; Based on the radiation field dose cutoff radius, multiple threads are invoked in parallel to determine the voxels of interest from the tissue voxels in the tissue voxel grid, including: For any one of the multiple tissue voxels, the distance between the tissue voxel and each of the multiple candidate voxels is calculated to obtain multiple distance values, wherein the candidate voxels represent tissue voxels traversed by the pencil beam; Determine the minimum distance value from the plurality of distance values; When the minimum distance value is less than the field dose cutoff radius, the tissue voxel and the candidate voxel corresponding to the minimum distance value are both determined as the voxel of interest, wherein the field dose cutoff radius characterizes the upper limit of the action distance of all the pencil beams in a field that make a significant contribution to the dose distribution in the lateral direction; Based on the dose kernel parameters of each of the multiple intersection points, multiple threads are invoked to execute the dose calculation tasks of multiple voxels of interest in parallel to obtain the voxel deposition dose of the voxels of interest. Based on the voxel deposition dose of each of the multiple voxels of interest, the dose distribution results associated with the specified tissue are determined.

2. The method according to claim 1, characterized in that, The pencil beam travels in the magnetic field; wherein, calculating the pencil beam trajectory of each of the multiple pencil beams includes: Multiple threads are used to compute the pencil beam trajectories of multiple pencil beams in parallel, wherein any one of the multiple threads is used to perform the following operation: When the velocity value of the pencil beam at the current intersection point is greater than a preset velocity threshold, the point magnetic field vector of the current intersection point where the pencil beam is located is determined from the three-dimensional magnetic field distribution map, and the three-dimensional magnetic field distribution map represents the spatial magnetic field generated by the magnetic resonance device. Based on the point magnetic field vector, the velocity vector of the pencil beam at the current intersection point, and the straight-line travel distance of the pencil beam within the voxel, the upper limit value of the deflection angle of the pencil beam at the current intersection point is determined, and the velocity vector includes the velocity value; The point magnetic field vector and the velocity vector are processed according to the intersection calculation algorithm that matches the upper limit of the deflection angle to obtain the next intersection point. The pencil beam trajectory includes multiple intersection points. The pencil-shaped beam trajectory is generated based on multiple intersection points.

3. The method according to claim 2, characterized in that, The intersection point calculation algorithm is determined based on the following operations: If the upper limit of the deflection angle is less than the preset deflection angle threshold, the intersection point calculation algorithm is determined to be the straight line intersection algorithm; If the upper limit of the deflection angle is greater than or equal to the preset deflection angle threshold, the intersection calculation algorithm is determined to be the spiral intersection algorithm.

4. The method according to claim 2, characterized in that, The velocity vector of the pencil beam at the current intersection point is determined based on the following operation: Based on the deflection angle, point magnetic field vector, and velocity vector of the pencil beam at the previous intersection point, the velocity vector of the pencil beam at the current intersection point is determined.

5. The method according to claim 1, characterized in that, The process involves calling multiple threads to execute dose calculation tasks for multiple voxels of interest in parallel, based on the dose kernel parameters of each of multiple intersection points, to obtain the voxel deposition dose of the voxels of interest, including: The following operation is performed using any one of the multiple threads for any given stroke bundle: Based on the foot of the perpendicular from the voxel center point of the voxel of interest to the trajectory of the pen beam, determine the radial distance from the voxel center point to the foot of the perpendicular and the path length of the pen beam from the exit point to the foot of the perpendicular. When the path length is less than the range of the pencil beam and the radial distance is less than the dose cutoff radius of the pencil beam, the path length and the dose kernel parameters of the pencil beam at multiple intersection points are processed using an interpolation algorithm to obtain the dose kernel parameters corresponding to the perpendicular foot. The dose of the pencil beam deposited on the voxel of interest is obtained based on the dose kernel parameters of the foot of the pendulum; The voxel deposition dose corresponding to the voxel of interest is obtained based on the dose of each of the multiple pencil beams deposited on the voxel of interest. Based on multiple threads, the voxel deposition dose of each of the multiple voxels of interest is obtained in parallel.

6. The method according to claim 5, characterized in that, The dose nuclear parameters include the Gaussian standard deviation of nuclear interactions, the Gaussian standard deviation of multiple Coulomb scattering, and the nuclear interaction weights; The dose kernel parameters of the pencil beam at each of the multiple intersection points are obtained based on the following operation: For any one of the plurality of intersection points, Based on the target beam size and the first beam size, the Gaussian standard deviation of the multiple Coulomb scattering is determined, wherein the target beam size characterizes the lateral broadening caused by the pencil beam traveling in the voxel, and the first beam size characterizes the lateral broadening of the pencil beam in the tissue voxel due to multiple Coulomb scattering. The first beam size is determined based on the path length of the pencil beam from the initial tissue intersection point to the current intersection point, and the initial tissue intersection point characterizes the first intersection point between the pencil beam and the tissue voxel mesh. The Gaussian standard deviation of the nuclear interaction is determined based on the Gaussian standard deviation of the multiple Coulomb scattering and the second beam spot size, wherein the second beam spot size characterizes the lateral broadening of the pencil beam in the tissue voxel due to nuclear interaction, and the second beam spot size is determined based on the path length of the pencil beam from the initial tissue intersection point to the current intersection point; The nuclear interaction weights are determined based on the path length of the pencil beam from the initial tissue intersection point to the current intersection point.

7. The method according to claim 6, characterized in that, The target beam size is obtained based on the following operation: Based on the initial energy of the pencil beam, the strength of the main magnetic field, and the direction of the magnetic field configuration, the initial optical parameters are looked up from the optical parameter lookup table. The target beam spot size is determined based on the initial optical parameters and the path length of the pencil beam from the exit point to the current intersection point.

8. The method according to claim 1, characterized in that, The field dose cutoff radius is obtained based on the following operation: For any one of the pencil beams in a single firing field The pencil beam trajectory is sampled based on a preset sampling rule to obtain multiple sampling points; Based on the lateral dose distribution of the sampling points, the dose cutoff radius corresponding to the sampling points is obtained by using the lateral dose integral function, wherein the lateral dose distribution is obtained by interpolating the dose kernel parameters of each of the multiple intersection points on the pencil beam trajectory; The maximum dose cutoff radius among the dose cutoff radii of each of the multiple sampling points is determined as the dose cutoff radius of the pencil beam; The maximum dose cutoff radius among the dose cutoff radii of the plurality of pencil beams is determined as the field dose cutoff radius.

9. A proton dose distribution determination device based on a graphics processor, comprising: The trajectory calculation module is used to calculate the trajectory of each of the multiple pencil beams, wherein the intersection points in the pencil beam trajectory include the intersection points where the pencil beam intersects with the medium voxel grid of the transmission medium, and the intersection points where the pencil beam intersects with the tissue voxel grid of a specified tissue, wherein the transmission medium is air; The acquisition module is used to acquire the dose kernel parameters of the pencil beam at multiple intersection points; The first determining module is used to call multiple threads in parallel to determine the voxel of interest from the tissue voxels in the tissue voxel grid according to the field dose cutoff radius. The field dose cutoff radius represents the upper limit of the action distance of all the pencil beams under a field that make a significant contribution to the dose distribution in the lateral direction. The first determining module includes a first determining submodule, a second determining submodule, and a third determining submodule. The first determination submodule is used to calculate the distance between any one of the multiple tissue voxels and each of the multiple candidate voxels to obtain multiple distance values, wherein the candidate voxels represent tissue voxels that are traversed by the pencil beam. The second determination submodule is used to determine the minimum distance value from multiple distance values; The third determination submodule is used to determine both the tissue voxel and the candidate voxel corresponding to the minimum distance value as the voxel of interest when the minimum distance value is less than the cutoff radius of the radiation field dose. The calling module is used to call multiple threads to execute dose calculation tasks of multiple voxels of interest in parallel based on the dose kernel parameters of each of the multiple intersection points, so as to obtain the voxel deposition dose of the voxels of interest; The second determining module is used to determine the dose distribution result related to a specified tissue based on the voxel deposition dose of each of the multiple voxels of interest.

Images