Adaptive ray-tracing based front projection method using gpu parallelization

Abstract

This invention provides a GPU-parallelized adaptive ray trajectory pre-projection method, comprising: determining the geometric parametric equations of the ray; each GPU thread executing single-ray calculations; for each single ray, determining whether the single ray passes through the valid range of the 3D volume data; if so, the single ray is valid; for valid single rays, adaptively selecting the main step axis based on the slope of the single ray, and filtering to obtain valid sampling points; calculating the voxel index and sub-voxel offset, and obtaining the linear attenuation coefficient through bilinear interpolation sampling; performing path length weighted summation on all valid sampling points to obtain the grayscale value of the projected pixel; and integrating to obtain the pre-projection result. This invention's method eliminates invalid rays through 3D boundary pre-judgment, adaptively selects the main step axis based on the ray slope, achieves high-precision integration by combining bilinear interpolation and path length correction, and achieves synchronous calculation through a parallel architecture, significantly improving computational efficiency.

Description

Technical Field

[0001] This invention belongs to the field of medical image reconstruction and GPU parallel computing technology. Specifically, it relates to an adaptive tracking ray trajectory forward projection method based on GPU parallelization. It is applicable to scenarios such as digital reconstructed radiographic images (DRR) generation, iterative reconstruction forward projection calculation, image registration and quality assurance of medical imaging equipment such as CBCT and CT. It can significantly improve the forward projection calculation speed while ensuring imaging accuracy, and provide technical support for clinical intraoperative navigation, online reconstruction and equipment calibration. Background Technology

[0002] Preprojection calculation is a core foundational algorithm in the field of medical imaging. It generates digitally reconstructed radiographic images (DRR) in the form of two-dimensional projected images by simulating the attenuation process of X-rays passing through three-dimensional volume data. It is widely used in scenarios such as CBCT image reconstruction, radiotherapy planning, 2D-3D image registration, and equipment geometric calibration. The accuracy and speed of preprojection calculation directly determine the quality of image reconstruction and the efficiency of clinical work.

[0003] Existing front projection technology mainly suffers from the following core problems:

[0004] 1. CPU serial computing efficiency is extremely low: Traditional CPU-based front projection algorithms use a ray-by-ray serial computing mode. When faced with a projection plane with millions of pixels and three-dimensional data with tens of millions of voxels, the calculation time for a single frame of projection can reach several seconds to tens of seconds, which cannot meet the speed requirements of scenarios such as intraoperative real-time navigation and online iterative reconstruction.

[0005] 2. Traditional ray tracing algorithms struggle to balance accuracy and speed: The classic Siddon algorithm uses the intersection of rays and voxel surfaces, which is highly accurate but the intersection calculation is complex and has many branches. In GPU parallel scenarios, it is prone to thread divergence and low parallel efficiency. The simplified equidistant stepping algorithm is fast, but the uneven distribution of sampling points significantly increases the error in the oblique ray direction, making it difficult to meet clinical requirements for imaging quality.

[0006] 3. Existing GPU parallel solutions are not optimized enough: Some GPU front projection methods rely on the GPU's built-in hardware texture interpolation for acceleration, which is limited by the hardware interpolation accuracy and memory bandwidth, and cannot flexibly customize the interpolation and correction logic; other algorithms only perform simple thread mapping and do not design adaptive calculation strategies for ray geometry characteristics, resulting in a large number of invalid calculations and redundant branches, and low GPU computing power utilization.

[0007] 4. Poor versatility and configurability: Most algorithms are designed for specific volume data sizes and fixed projection geometry, and cannot flexibly adapt to phantom / human body data with different fields of view, CBCT systems with different source axis distances and source image distances, and it is also difficult to support batch calculation of multi-angle sequence projections.

[0008] Therefore, there is an urgent need for a fast forward projection method that balances high precision and high speed, is based on a GPU fully parallel architecture, and has an adaptive ray tracing mechanism, to solve the core pain points of existing technologies, such as long computation time, low parallel efficiency, and the inability to achieve both precision and speed. Summary of the Invention

[0009] This invention aims to provide a GPU-parallelized adaptive ray trajectory tracking front projection method to solve the problems of slow computation speed, low GPU parallel efficiency, large oblique ray sampling error, and poor versatility of existing front projection algorithms.

[0010] To achieve the above objectives, this invention provides a GPU-parallelized adaptive ray trajectory forward projection method, comprising:

[0011] Step S1: Initialize the system's constant parameters and configure GPU resources;

[0012] Step S2: Preprocess and store the 3D volume data;

[0013] Step S3: Generate the coordinate array of the projection plane and determine the geometric parametric equations of the ray;

[0014] Step S4: GPU parallelization is adopted, and each parallel GPU thread executes the calculation of a single ray. For each single ray, the coordinates of its intersection point with the boundary plane of the 3D volume data are pre-calculated to determine whether the single ray passes through the valid range of the 3D volume data. If so, the single ray is valid.

[0015] Step S5: For a valid single ray, adaptively select the main stepping axis based on the slope of the single ray, and step-track the single ray to obtain the coordinates of the sampling point outside the main stepping axis, and filter out the valid sampling points.

[0016] Step S6: For each valid sampling point, calculate the voxel index and sub-voxel offset of the valid sampling point in the 3D volume data, and obtain the linear attenuation coefficient of the sampling point location by bilinear interpolation sampling; then, perform path length weighted summation on all valid sampling points on a single ray to obtain the gray value of the projection pixel corresponding to the single ray; integrate the gray values ​​of the projection pixels corresponding to all valid single rays to obtain the front projection result.

[0017] The constant parameters of the initialization system are determined based on the input imaging geometry parameters and the parameters of the 3D volume data; the constant parameters of the system include: the number of rows of the projected image. Number of columns of projected images Dimensions of 3D volume data in the X / Y / Z directions Voxel spacing of 3D volume data Distance from light source to rotation center (SAD), distance from light source to detector (SID); total number of projection angles and angle intervals;

[0018] Configure GPU resources, specifically including: writing the system's constant parameters into the GPU's constant memory, allocating volume data storage space, projected image storage space, and coordinate array storage space in the GPU's global memory, and completing the GPU memory initialization and zeroing.

[0019] In step S2, the CT values ​​of the three-dimensional volume data are converted into X-ray linear attenuation coefficients through linear transformation to preprocess the three-dimensional volume data.

[0020] In step S3, determining the geometric parameter equation of the ray specifically includes: for each projection angle, based on the initial position of the light source and the initial position of the ray endpoint, calculating the transformed ray start point and the transformed ray endpoint after rotating the projection angle, and determining the geometric parameter equation of a single ray.

[0021] In step S4, the effective physical range of the three-dimensional volume data is:

[0022] ,

[0023] The formula for calculating boundary extrema is:

[0024] ,

[0025] Where x, y, and z are the x, y, and z coordinates of the three-dimensional volume data, respectively. , , These are the boundary extrema of the three-dimensional volume data in the x, y, and z directions, respectively. , , These are the dimensions of the 3D volume data in the x, y, and z directions, respectively. These are the voxel spacings of the three-dimensional volume data in the x, y, and z directions, respectively.

[0026] By pre-calculating the coordinates of the intersection point between the ray and the boundary plane of the 3D volume data, it is determined whether the single ray passes through the valid range of the 3D volume data. Specifically, this includes: obtaining the coordinates of the intersection point between the single ray and the boundary plane of the 3D volume data. , , If the coordinates of the intersection points of the planes with the six boundary extrema fall within the valid range of the corresponding coordinate axes, then the single ray is determined to intersect the 3D volume data, and the subsequent step tracking process can be initiated; otherwise, the projection value of the single ray is recorded as 0, and the GPU thread calculation corresponding to the single ray is terminated directly.

[0027] Step S5 specifically includes:

[0028] Step S51: For a valid single ray, calculate the slope of the single ray in the XOZ plane based on the two coordinate values ​​of the single ray. ;

[0029] Step S52: Adaptively select the main stepping axis based on the absolute value of the slope of the single ray in the XOZ plane to ensure that the sampling points on the single ray are evenly distributed along the ray direction;

[0030] Step S53: Make the sampling points move at equal intervals along the main stepping axis; each step is one step size, and calculate the coordinate values ​​of the current sampling point on the two coordinate axes other than the main stepping axis according to the ray parameter equation;

[0031] Step S54: For each sampling point, determine whether the coordinate values ​​of the sampling point on the two coordinate axes other than the main step axis fall within the valid range of the three-dimensional volume data. If so, determine it as a valid sampling point and proceed to step S6.

[0032] In step S52, if the slope of the single ray in the XOZ plane Choose the Z-axis as the main stepping axis, so that the sampling points move along the Z-axis from... arrive With step size Equal interval stepping, The boundary extrema of the 3D volume data in the z-direction;

[0033] If the slope of the single ray in the XOZ plane Choose the X-axis as the primary stepping axis, so that the sampling points move along the X-axis from... arrive With step size Equal interval stepping, represents the boundary extremum in the x-direction of the 3D volume data.

[0034] In step S6, when the Z-axis is the main stepping axis, the voxel index and sub-voxel offset weights of the effective sampling points in the X and Y directions are:

[0035] ,

[0036] ,

[0037] in, The voxel index of the effective sampling point in the X direction. The voxel index of the effective sampling point in the Y direction. The subvoxel offset weights of the effective sampling points in the X direction, The subvollar offset weights for effective sampling points in the Y direction, where x is the x-coordinate of the sampling point and y is the y-coordinate of the sampling point. , These are the boundary extrema of the three-dimensional volume data in the x and y directions, respectively. These are the voxel spacings in the x and y directions of the three-dimensional volume data, respectively.

[0038] Linear attenuation coefficient at the location of the sampling point for:

[0039] ,

[0040] in, , , , For four adjacent voxel values, , , These are the voxel indices of the valid sampling points in the X, Y, and Z directions, respectively. The subvoxel offset weights of the effective sampling points in the X direction, The subvoxel offset weights of the effective sampling points in the Y direction.

[0041] The grayscale value of the projected pixel corresponding to the single ray is obtained by weighting and summing the path lengths of all valid sampling points on the single ray, specifically including:

[0042] Step S61: Introduce a dual path length correction factor to correct the spindle step size to the actual ray path length. ;

[0043] The dual path length correction factor includes: XOZ plane tilt correction factor and three-dimensional direction correction coefficient ;

[0044] After spindle step size correction, the actual ray path length for:

[0045] ,

[0046] in, For three-dimensional orientation correction coefficients; Let be the slope of the single ray in the XOZ plane. It is the voxel spacing in the z-direction of the 3D volume data;

[0047] Three-dimensional orientation correction coefficient The calculation formula is:

[0048] ,

[0049] in, These are the x, y, and z coordinates of the ray's origin, respectively. These are the x, y, and z coordinates of the ray's endpoint, respectively.

[0050] Step S62: Based on the actual path length of the ray The grayscale value of the projected pixel corresponding to the single ray is obtained by weighting and summing the path lengths of all valid sampling points on the single ray.

[0051] The GPU-parallelized adaptive ray trajectory pre-projection method further includes step S7: traversing all projection angles and repeating steps S3 to S6 to obtain the pre-projection results for all projection angles.

[0052] The method of this invention adopts a ray-driven calculation mode, eliminates invalid rays by pre-judging the three-dimensional boundary, adaptively selects the stepping principal axis based on the ray slope, achieves high-precision integration by combining bilinear interpolation and path length correction, and achieves synchronous calculation of massive rays by relying on the GPU parallel architecture. It significantly improves the calculation efficiency while ensuring projection accuracy and is compatible with the front projection and DRR generation requirements of various CBCT / CT equipment. Attached Figure Description

[0053] Figure 1 This is a flowchart of a method for providing an adaptive ray trajectory tracking forward projection based on GPU parallelization according to the present invention;

[0054] Figures 2A-2C These are the three views of the input 3D volume data. The 3D volume data can be either CT data or CBCT data. Figures 2A-2C It shows three different perspectives;

[0055] Figure 3A It is a front projection image obtained by applying the front projection method of this patent. Figure 3B These are actual X-ray images taken in DR mode;

[0056] Figure 4A and Figure 4B This is a comparison diagram of the effects of the traditional front projection method and the front projection method of this patent, in which... Figure 4A This is a topographic image of the spherical marker points inside the head mold under the traditional front projection method. Figure 4B This is a topographic diagram of the spherical marker points inside the head mold under the front projection method of the present invention. Detailed Implementation

[0057] The technical principle of the GPU-parallelized adaptive ray trajectory tracking pre-projection method of this invention is a fully parallel processing architecture: "parameter initialization → volume data preprocessing → ray geometry calculation → boundary intersection determination → adaptive step tracking → bilinear interpolation integration → multi-view projection output". Relying on the CUDA parallel computing model, each projected pixel corresponds to a single independent computation thread, achieving complete parallelization of ray trajectory tracking and attenuation integration.

[0058] like Figure 1 As shown, based on the above principles, the present invention provides a GPU-parallelized adaptive ray trajectory pre-projection method that specifically includes the following steps:

[0059] Step S1: Initialize the system's constant parameters and configure GPU resources;

[0060] The constant parameters for initializing the system are determined based on the input imaging geometry parameters and the parameters of the 3D volume data. These constant parameters include: the number of rows in the projected image. Number of columns of projected images Dimensions of 3D volume data in the X / Y / Z directions Voxel spacing of 3D volume data ; Distance from light source to rotation center SAD, distance from light source to detector SID; Total number of projection angles and angle interval.

[0061] Configure GPU resources, specifically including: writing the system's constant parameters into the GPU's constant memory, allocating volume data storage space, projected image storage space, and coordinate array storage space in the GPU's global memory, and completing the GPU memory initialization and zeroing.

[0062] The GPU's memory includes constant memory shared by all GPU threads, global memory that can be read and written by all GPU threads, registers and local memory exclusive to each GPU thread. In this embodiment, constant memory is used to store the system's constant parameters; global memory is used to store 3D volume data, the pre-projection results of each frame (i.e., digital reconstructed radiographic images DRR), and coordinate arrays; registers are used to store the x / z coordinates, slope, interpolation weights, and accumulated projection values ​​of each single ray; and local memory is used to store temporary variables that cannot be stored in the registers.

[0063] Step S2: Preprocess and store the 3D volume data;

[0064] The CT values ​​of the 3D volume data input from the host are transmitted to the GPU device, which is usually the CPU, to obtain the 3D volume data.

[0065] Figures 2A-2CThese are the three views of the input 3D volume data. The 3D volume data can be either CT data or CBCT data. Figures 2A-2C It shows three different perspectives.

[0066] The CT values ​​of 3D volume data are converted into X-ray linear attenuation coefficients through linear transformation to preprocess the 3D volume data. The conversion formula for CT values ​​to X-ray linear attenuation coefficients is as follows:

[0067] ,

[0068] in, These are the original CT values. This is the scaling factor. This is the linear attenuation coefficient after conversion.

[0069] The preprocessed 3D volume data is stored sequentially in the GPU's global memory in the ZYX order for subsequent ray sampling.

[0070] Step S3: Generate the coordinate array of the projection plane and determine the geometric parametric equations of the ray;

[0071] Generate the coordinate array of the projection plane, with the origin located at the center of the projection plane. The formula for calculating the coordinate array of the projection plane is:

[0072] ,

[0073] ,

[0074] Where U represents the horizontal direction of the detector, corresponding to a complete row of detectors, and V represents the vertical direction of the detector, corresponding to a complete column of detectors. , where are the pixel spacing in the U-direction and V-direction of the detector, respectively, and j is the column index of the detector. Let be the column coordinate of the detector, and i be the row index of the detector. Let C be the row coordinate of the detector, C be the total number of columns of the detector, and R be the total number of rows of the detector.

[0075] Determining the geometric parametric equations of a ray involves, for each projection angle, calculating the transformed ray start point (source position) and transformed ray end point after rotating the projection angle, based on the initial position of the light source and the initial position of the ray endpoint, and then determining the geometric parametric equations of the individual ray. The initial position is the position when the projection angle is 0.

[0076] Among them, for the first Each projection angle Based on the initial position of the light source By rotating the first axis around the Y-axis Each projection angle The transformation yields the starting point of the transformed ray. , These are the x, y, and z coordinates of the ray's origin; and the pixels on the projection plane. The initial position of the corresponding ray endpoint is located on the detector plane. Based on the initial position of the ray endpoint, the transformed ray endpoint is obtained after rotation transformation. , These are the x, y, and z coordinates of the ray's endpoint, respectively. The distance from the light source to the detector. Let be the distance from the light source to the center of rotation. Finally, the geometric parametric equation for a single ray is expressed as:

[0077] ,

[0078] in, The starting point of the ray after the transformation, The endpoint of the transformed ray is represented by t, which is a time coefficient. 0 indicates the starting time of the ray, and 1 indicates the ending time of the current ray.

[0079] Step S4: GPU parallelization is adopted, and each parallel GPU thread executes the calculation of a single ray. For each single ray, the coordinates of its intersection point with the boundary plane of the 3D volume data are pre-calculated to determine whether the single ray passes through the valid range of the 3D volume data. If so, the single ray is valid.

[0080] This enables pre-determination of the boundary intersection between a single ray and three-dimensional volume data, eliminating invalid rays that do not intersect at all, and avoiding subsequent invalid calculations.

[0081] The effective physical range of the three-dimensional volume data is as follows:

[0082] ,

[0083] The formula for calculating boundary extrema is:

[0084] ,

[0085] Where x, y, and z are the x, y, and z coordinates of the three-dimensional volume data, respectively. , , These are the boundary extrema of the three-dimensional volume data in the x, y, and z directions, respectively. , , These are the dimensions of the 3D volume data in the x, y, and z directions, respectively. These are the voxel spacings of the three-dimensional volume data in the x, y, and z directions, respectively.

[0086] By pre-calculating the coordinates of the intersection point between the ray and the boundary plane of the 3D volume data, it is determined whether the single ray passes through the valid range of the 3D volume data. Specifically, this includes: obtaining the coordinates of the intersection point between the single ray and the boundary plane of the 3D volume data. , , If the coordinates of the intersection points (including three coordinate values) of the planes of the six boundary extrema fall within the valid range of the corresponding coordinate axes, then the single ray is determined to intersect the 3D volume data, and the subsequent step tracking process can be entered; otherwise, the projection value of the single ray is recorded as 0, and the GPU thread calculation corresponding to the single ray is directly terminated.

[0087] Step S5: For a valid single ray, adaptively select the main stepping axis based on the slope of the single ray, and step-track the single ray to obtain the coordinates of the sampling point outside the main stepping axis, and filter out the valid sampling points.

[0088] Step S5 specifically includes:

[0089] Step S51: For a valid single ray, calculate the slope of the single ray in the XOZ plane based on the two coordinate values ​​of the single ray. ;

[0090] Calculate the slope of a single ray in the XOZ plane. The reason is that, in a fixed coordinate system, the Z-axis connects the X-ray source to the detector plate. When the frame containing the X-ray source and detector plate rotates, the plane formed by a single ray is essentially the XOZ plane. Using this plane to calculate the slope is more reasonable and has better adaptability to various scenarios. The Y-axis can be understood as an extension axis from two dimensions to three dimensions, and is generally not considered.

[0091] The slope of a single ray in the XOZ plane for:

[0092] ,

[0093] in, Let z be the z-coordinate of the ray's origin. Let z be the z-coordinate of the ray's endpoint. Let x be the x-coordinate of the ray's origin. Let x be the x-coordinate of the ray's endpoint.

[0094] Step S52: Adaptively select the main stepping axis based on the absolute value of the slope of the single ray in the XOZ plane to ensure that the sampling points on the single ray are evenly distributed along the ray direction and that there is no missed sampling distortion.

[0095] Wherein, if the slope of a single ray in the XOZ plane Choose the Z-axis as the primary stepping axis (its interval is relatively smaller compared to choosing the X-axis), so that the sampling points move along the Z-axis from... arrive With step size Equal interval stepping, The boundary extrema of the 3D volume data in the z-direction;

[0096] If the slope of the single ray in the XOZ plane Choose the X-axis as the primary stepping axis (its interval is relatively smaller compared to choosing the Z-axis), so that the sampling points move along the X-axis from... arrive With step size Equal interval stepping, represents the boundary extremum in the x-direction of the 3D volume data.

[0097] Step S53: Make the sampling points move at equal intervals along the main stepping axis; each step is one step size, and calculate the coordinate values ​​of the current sampling point on the two coordinate axes other than the main stepping axis according to the ray parameter equation.

[0098] For example, when the Z-axis is the main stepping axis, the sampling points The x-coordinate corresponding to the coordinate y coordinate for:

[0099] ,

[0100] ,

[0101] in, Let z be the z-coordinate of the ray's origin. Let z be the z-coordinate of the ray's endpoint. Let x be the x-coordinate of the ray's origin. Let x be the x-coordinate of the ray endpoint and z be the z-coordinate of the sampling point. coordinate, Here is the x-coordinate of the sampling point. Here is the y-coordinate of the sampling point.

[0102] Step S54: For each sampling point, determine whether the coordinate values ​​of the sampling point on the two coordinate axes other than the main step axis fall within the valid range of the three-dimensional volume data. If so, determine it as a valid sampling point and proceed to step S6.

[0103] In this embodiment, if the x and y coordinates of a certain sampling point both fall within the valid range of the three-dimensional volume data, then the interpolation sampling and integration accumulation stage is entered.

[0104] Step S6: For each valid sampling point, calculate the voxel index and sub-voxel offset of the valid sampling point in the 3D volume data, and obtain the linear attenuation coefficient of the sampling point location by bilinear interpolation sampling; then, perform path length weighted summation on all valid sampling points on a single ray to obtain the gray value of the projection pixel corresponding to the single ray; integrate the gray values ​​of the projection pixels corresponding to all valid single rays to obtain the front projection result.

[0105] When the Z-axis is the main stepping axis, the voxel index and sub-voxel offset weights of the effective sampling points in the X and Y directions are:

[0106] ,

[0107] ,

[0108] in, The voxel index of the effective sampling point in the X direction. The voxel index of the effective sampling point in the Y direction. The subvoxel offset weights of the effective sampling points in the X direction, The subvollar offset weights for effective sampling points in the Y direction, where x is the x-coordinate of the sampling point and y is the y-coordinate of the sampling point. , These are the boundary extrema of the three-dimensional volume data in the x and y directions, respectively. These are the voxel spacings in the x and y directions of the three-dimensional volume data, respectively.

[0109] Voxel index of valid sampling points in the Z direction for:

[0110] ,

[0111] Its step length The steps are spaced evenly, so rounding is not required.

[0112] The linear attenuation coefficient at the location of the sampling point is obtained by bilinear interpolation sampling. Specifically, it includes: taking four adjacent voxel values ​​and performing bilinear interpolation based on the voxel index and sub-voxel offset of the effective sampling point in the three-dimensional volume data to obtain the linear attenuation coefficient at the location of the sampling point.

[0113] Among them, the linear attenuation coefficient at the location of the sampling point for:

[0114] ,

[0115] in, , , , For four adjacent voxel values, , , These are the voxel indices of the valid sampling points in the X, Y, and Z directions, respectively. The subvoxel offset weights of the effective sampling points in the X direction, The subvoxel offset weights of the effective sampling points in the Y direction.

[0116] Specifically, the grayscale value of the projected pixel corresponding to the single ray is obtained by weighting and summing the path lengths of all valid sampling points on the single ray, including:

[0117] Step S61: Introduce a dual path length correction factor to correct the spindle step size to the actual ray path length. ;

[0118] The dual path length correction factor includes: XOZ plane tilt correction factor. and three-dimensional direction correction coefficient .

[0119] Therefore, after the spindle step size is corrected, the actual ray path length is... for:

[0120] ,

[0121] in, The three-dimensional orientation correction coefficient characterizes the proportional relationship between the projection of the ray onto the XOZ plane and the total length of the ray. Let be the slope of the single ray in the XOZ plane. It is the voxel spacing in the z-direction of the 3D volume data, i.e., the principal axis step size.

[0122] Three-dimensional orientation correction coefficient The calculation formula is:

[0123] ,

[0124] in, These are the x, y, and z coordinates of the ray's origin, respectively. These are the x, y, and z coordinates of the ray's endpoint, respectively.

[0125] Step S62: Based on the actual path length of the ray The grayscale value of the projected pixel corresponding to the single ray is obtained by weighting and summing the path lengths of all valid sampling points on the single ray.

[0126] Wherein, the contribution of each valid sampling point to the projected value is Therefore, by performing a path length-weighted summation on all valid sampling points along a single ray, the grayscale value of the projected pixel corresponding to the single ray is obtained. for:

[0127] ,

[0128] in, The linear attenuation coefficient at the location of the sampling point. This represents the actual path length of the ray.

[0129] Step S6 can also employ simplified stepping and interpolation logic for specific scenarios to further improve computational efficiency. For example: 1) When the ray is perpendicular to the coordinate plane, i.e., the ray path is along the Z-axis, and only the Z-value changes during the penetration process, then only pixel accumulation along the Z-axis is needed. 2) When the ray penetration path only changes in two dimensions, such as X and Z or Y and Z, then three-dimensional interpolation and correction can be reduced to two-dimensional interpolation and correction, reducing redundant computation. Step S7: Traverse all projection angles and repeat steps S3 to S6 to obtain the pre-projection results for all projection angles.

[0130] Each time a projection angle is completed, the projection data from the GPU is copied back to the host for storage, and the GPU's projection buffer is reset. After all projection angles have been calculated, the GPU's video memory resources are released, and the pre-projection results of all projection angles are output as a complete multi-angle projection sequence.

[0131] In this embodiment, the values ​​and functions of the core parameters are shown in Table 1. In other embodiments, these values ​​can be changed as needed.

[0132] Table 1: Explanation of Core Parameters

[0133] Experimental results:

[0134] Step S1: Initialize the system's constant parameters and configure GPU resources;

[0135] Input configuration parameters: Projected image 1536×1536 pixels, volume data 512×512×512 voxels, voxel pitch 0.5mm, SAD=2100mm, SID=3000mm, projection angles 360, scale factor 0.0171. Call the initialization interface to set up the GPU device, allocate video memory, and write geometric constants to the GPU's constant memory.

[0136] Step S2: Preprocess and store the 3D volume data;

[0137] The 512×512×512 CT value volume data from the host machine is copied to the GPU memory. A 3D parallel kernel function is started to perform the conversion of CT values ​​to linear attenuation coefficients on a voxel-by-voxel basis, relying on the formula... Output the preprocessed attenuation coefficient volume data.

[0138] Step S3: Generate the coordinate array of the projection plane and determine the geometric parametric equations of the ray;

[0139] On the host side, based on the calculation formulas of U and V coordinates, the physical coordinate array is calculated and copied to the GPU memory; the initial position of the light source and the initial position of the ray endpoint are calculated, and the coordinates of the rotated light source and the ray endpoint are calculated based on the current projection angle, thereby determining the geometric parameter equation of the ray.

[0140] Step S4: GPU parallelization is adopted, and each parallel GPU thread executes the calculation of a single ray. For each single ray, the coordinates of its intersection point with the boundary plane of the 3D volume data are pre-calculated to determine whether the single ray passes through the valid range of the 3D volume data. If so, the single ray is valid.

[0141] Step S5: For a valid single ray, adaptively select the main stepping axis based on the slope of the single ray, and step-track the single ray to obtain the coordinates of the sampling point outside the main stepping axis, and filter out the valid sampling points.

[0142] Step S6: For each valid sampling point, calculate the voxel index and sub-voxel offset of the valid sampling point in the 3D volume data, and obtain the linear attenuation coefficient of the sampling point location by bilinear interpolation sampling; then, perform path length weighted summation on all valid sampling points on a single ray to obtain the gray value of the projection pixel corresponding to the single ray; integrate the gray values ​​of the projection pixels corresponding to all valid single rays to obtain the front projection result.

[0143] In other words, in steps S4-S6, a two-dimensional parallel kernel function is started, with each thread corresponding to a single ray and thus a projected pixel, and ray tracing is performed independently: first, a three-dimensional boundary intersection determination is performed, and invalid rays are directly assigned a value of 0; for valid rays, the XOZ plane slope is calculated. The stepping axis is adaptively selected, and the stepping is performed at equal intervals along the main axis. Bilinear interpolation and path length weighted accumulation are performed point by point.

[0144] Figure 3A It is a front projection image obtained by applying the front projection method of this patent. Figure 3B These are actual X-ray images taken in DR mode; comparing the two images from the same angle, it can be seen that the result of the front projection is basically consistent with the actual image.

[0145] Figure 4A and Figure 4BThis is a comparison diagram of the effects of the traditional front projection method and the front projection method of this patent, in which... Figure 4A This is a topographic image of the spherical marker points inside the head mold under the traditional front projection method. Figure 4B This is a topographical diagram of the spherical marker points inside the head mold under the front projection method of the present invention. From Figure 4A and Figure 4B It can be seen that the front projection method of the present invention processes the image more smoothly and naturally, and also conforms to the objective reality of the sphere's shape.

[0146] Step S7: Traverse all projection angles and repeat steps S3 to S6 to obtain the front projection results of all projection angles.

[0147] After the kernel function executes, the projection results are copied from the GPU back to the host and stored in the projection sequence, and the projection buffer is reset. After all angle calculations are completed, 360 frames of projection sequence data are output, GPU memory resources are released, and the entire front projection process is completed.

[0148] The adaptive ray trajectory tracking pre-projection method based on GPU parallelization of the present invention has the following advantages:

[0149] 1. Significantly improved computation speed: The method of this invention adopts a CUDA pixel-level fully parallel architecture, with one GPU thread corresponding to one ray. All projected pixels are calculated independently and synchronously, which can be more than 100 times faster than the CPU serial algorithm. The calculation of a single frame of high-definition projection can be completed in milliseconds (depending on the performance of the workstation, it takes 10-20ms per frame under normal working equipment, which meets the needs of engineering), and meets the requirements of real-time DRR generation and online iterative reconstruction.

[0150] 2. High and stable projection accuracy: The method of this invention adopts an adaptive principal axis stepping strategy, dynamically selecting the X or Z axis as the stepping direction based on the ray slope, avoiding the sparse sampling problem of oblique rays caused by single-direction stepping, and improving sampling uniformity while ensuring computational simplicity. In addition, a dual path length correction mechanism of "trigonometric function correction + three-dimensional scale correction" is adopted, which accurately corrects the principal axis step length to the actual path length of the ray based on the combined formula, ensuring the consistency of the integral of rays with different angles and slopes. Furthermore, combined with bilinear interpolation and the dual path length correction mechanism, the projection value integration is accurate, and the imaging quality of edges and oblique structures is excellent, with accuracy superior to traditional equal-interval stepping algorithms.

[0151] This invention achieves the accuracy of the near-surface intersection algorithm with extremely low computational overhead through a combination of adaptive stepping, bilinear interpolation, and dual path correction, while maintaining the high parallel efficiency of the equidistant stepping algorithm.

[0152] 3. High parallel efficiency and excellent computing power utilization: The method of this invention eliminates invalid rays in advance by pre-judging the three-dimensional boundary, reducing redundant calculations and forming a low-branch, high-parallelism front projection computing architecture; the core calculation process has simple branches, effectively reducing GPU thread divergence; constant memory stores geometric parameters and continuous memory layout data, giving full play to the GPU memory bandwidth advantage, and the computing power utilization is higher than that of conventional GPU front projection schemes.

[0153] 4. High versatility and adaptability: All core parameters (volume data size, resolution, imaging geometry, number of projection angles, etc.) can be flexibly configured, adapting to CBCT and CT equipment with different fields of view, supporting both single-frame DRR and multi-angle sequence projection modes, without the need to redesign the algorithm logic for specific hardware.

[0154] 5. Simplified Engineering Implementation: The method of this invention has a clear logic, no complex dependencies, and can be directly integrated into existing medical image reconstruction systems; it supports standard CUDA architecture GPUs, requiring no dedicated hardware; volume data size, resolution, and imaging geometry parameters can be flexibly adjusted according to the device model, supporting the front projection requirements of various devices such as dental CBCT, orthopedic CBCT, and whole-body CT; it also supports both single-frame and sequence working modes, adapting to different application scenarios such as calibration, navigation, and reconstruction, facilitating maintenance and secondary development. Compared to traditional CPU front projection algorithms, this method is more than two orders of magnitude faster; compared to conventional GPU front projection schemes, this method offers significantly better oblique structure imaging quality, better meeting the actual needs of real-time clinical DRR and online reconstruction, with engineering value far exceeding existing technologies.

[0155] The above description is merely a preferred embodiment of the present invention and is not intended to limit the scope of the invention. Various variations can be made to the above embodiments of the present invention. All simple and equivalent changes and modifications made in accordance with the claims and description of this application fall within the protection scope of the claims of this patent. All aspects not described in detail in this invention are conventional technical content.

Claims

1. A GPU-parallelized adaptive ray trajectory pre-projection method, characterized in that, include: Step S1: Initialize the system's constant parameters and configure GPU resources; Step S2: Preprocess and store the 3D volume data; Step S3: Generate the coordinate array of the projection plane and determine the geometric parametric equations of the ray; Step S4: GPU parallelization is adopted, and each parallel GPU thread executes the calculation of a single ray. For each single ray, the coordinates of its intersection point with the boundary plane of the 3D volume data are pre-calculated to determine whether the single ray passes through the valid range of the 3D volume data. If so, the single ray is valid. Step S5: For a valid single ray, adaptively select the main stepping axis based on the slope of the single ray, and step-track the single ray to obtain the coordinates of the sampling point outside the main stepping axis, and filter to obtain the valid sampling points; Step S6: For each valid sampling point, calculate the voxel index and sub-voxel offset of the valid sampling point in the 3D volume data, and obtain the linear attenuation coefficient of the sampling point location by bilinear interpolation sampling. Subsequently, the path length of all valid sampling points on a single ray is weighted and accumulated to obtain the grayscale value of the projected pixel corresponding to the single ray; the grayscale values ​​of the projected pixels corresponding to all valid single rays are integrated to obtain the front projection result.

2. The GPU-parallelized adaptive ray trajectory forward projection method according to claim 1, characterized in that, The constant parameters of the initialization system are determined based on the input imaging geometry parameters and the parameters of the 3D volume data; the constant parameters of the system include: the number of rows of the projected image. Number of columns of projected images Dimensions of 3D volume data in the X / Y / Z directions Voxel spacing of 3D volume data Distance from light source to rotation center (SAD), distance from light source to detector (SID); total number of projection angles and angle intervals; Configure GPU resources, specifically including: writing the system's constant parameters into the GPU's constant memory, allocating volume data storage space, projected image storage space, and coordinate array storage space in the GPU's global memory, and completing the GPU memory initialization and zeroing.

3. The GPU-parallelized adaptive ray trajectory forward projection method according to claim 1, characterized in that, In step S2, the CT values ​​of the three-dimensional volume data are converted into X-ray linear attenuation coefficients through linear transformation to preprocess the three-dimensional volume data.

4. The GPU-parallelized adaptive ray trajectory forward projection method according to claim 1, characterized in that, In step S3, determining the geometric parameter equation of the ray specifically includes: for each projection angle, based on the initial position of the light source and the initial position of the ray endpoint, calculating the transformed ray start point and the transformed ray endpoint after rotating the projection angle, and determining the geometric parameter equation of a single ray.

5. The GPU-parallelized adaptive ray trajectory forward projection method according to claim 1, characterized in that, In step S4, the effective physical range of the three-dimensional volume data is: ， The formula for calculating boundary extrema is: ， Where x, y, and z are the x, y, and z coordinates of the three-dimensional volume data, respectively. , , These are the boundary extrema of the three-dimensional volume data in the x, y, and z directions, respectively. , , These are the dimensions of the 3D volume data in the x, y, and z directions, respectively. These are the voxel spacings of the three-dimensional volume data in the x, y, and z directions, respectively. By pre-calculating the coordinates of the intersection point between the ray and the boundary plane of the 3D volume data, it is determined whether the single ray passes through the valid range of the 3D volume data. Specifically, this includes: obtaining the coordinates of the intersection point between the single ray and the boundary plane of the 3D volume data. , , If the coordinates of the intersection points of the planes with the six boundary extrema fall within the valid range of the corresponding coordinate axes, then the single ray is determined to intersect the 3D volume data, and the subsequent step tracking process can be initiated; otherwise, the projection value of the single ray is recorded as 0, and the GPU thread calculation corresponding to the single ray is terminated directly.

6. The GPU-parallelized adaptive ray trajectory forward projection method according to claim 1, characterized in that, Step S5 specifically includes: Step S51: For a valid single ray, calculate the slope of the single ray in the XOZ plane based on the two coordinate values ​​of the single ray. ; Step S52: Adaptively select the main stepping axis based on the absolute value of the slope of the single ray in the XOZ plane to ensure that the sampling points on the single ray are evenly distributed along the ray direction; Step S53: Make the sampling points move at equal intervals along the main stepping axis; each step is one step size, and calculate the coordinate values ​​of the current sampling point on the two coordinate axes other than the main stepping axis according to the ray parameter equation; Step S54: For each sampling point, determine whether the coordinate values ​​of the sampling point on the two coordinate axes other than the main step axis fall within the valid range of the three-dimensional volume data. If so, determine it as a valid sampling point and proceed to step S6.

7. The GPU-parallelized adaptive ray trajectory forward projection method according to claim 6, characterized in that, In step S52, if the slope of the single ray in the XOZ plane Choose the Z-axis as the main stepping axis, so that the sampling points move along the Z-axis from... arrive With step size Equal interval stepping, The boundary extrema of the 3D volume data in the z-direction; If the slope of the single ray in the XOZ plane Choose the X-axis as the primary stepping axis, so that the sampling points move along the X-axis from... arrive With step size Equal interval stepping, represents the boundary extremum in the x-direction of the 3D volume data.

8. The GPU-parallelized adaptive ray trajectory forward projection method according to claim 1, characterized in that, In step S6, when the Z-axis is the main stepping axis, the voxel index and sub-voxel offset weights of the effective sampling points in the X and Y directions are: ， ， in, The voxel index of the effective sampling point in the X direction. The voxel index of the effective sampling point in the Y direction. The subvoxel offset weights of the effective sampling points in the X direction, The subvollar offset weights for effective sampling points in the Y direction, where x is the x-coordinate of the sampling point and y is the y-coordinate of the sampling point. , These are the boundary extrema of the three-dimensional volume data in the x and y directions, respectively. These are the voxel spacings in the x and y directions of the three-dimensional volume data, respectively. Linear attenuation coefficient at the location of the sampling point for: ， in, , , , For four adjacent voxel values, , , These are the voxel indices of the valid sampling points in the X, Y, and Z directions, respectively. The subvoxel offset weights of the effective sampling points in the X direction, The subvoxel offset weights of the effective sampling points in the Y direction.

9. The GPU-parallelized adaptive ray trajectory forward projection method according to claim 1, characterized in that, The grayscale value of the projected pixel corresponding to the single ray is obtained by weighting and summing the path lengths of all valid sampling points on the single ray, specifically including: Step S61: Introduce a dual path length correction factor to correct the spindle step size to the actual ray path length. ; The dual path length correction factor includes: XOZ plane tilt correction factor and three-dimensional direction correction coefficient ; After spindle step size correction, the actual ray path length for: ， in, For three-dimensional orientation correction coefficients; Let be the slope of the single ray in the XOZ plane. It is the voxel spacing in the z-direction of the 3D volume data; Three-dimensional orientation correction coefficient The calculation formula is: ， in, These are the x, y, and z coordinates of the ray's origin, respectively. These are the x, y, and z coordinates of the ray's endpoint, respectively. Step S62: Based on the actual path length of the ray The grayscale value of the projected pixel corresponding to the single ray is obtained by weighting and summing the path lengths of all valid sampling points on the single ray.

10. The GPU-parallelized adaptive ray trajectory forward projection method according to claim 1, characterized in that, It also includes step S7: traversing all projection angles and repeating steps S3 to S6 to obtain the front projection results of all projection angles.

Images