A 2D / 3D image registration method and system based on Radon domain transformation

Abstract

The application discloses a 2D / 3D image registration method and system based on Radon domain transformation, and the registration method comprises two key components: a pre-processing technology based on Radon domain transformation and a neural network model prediction. The application adds 3D structure information in the network input, and trains the network, so that the problem of poor universality of the existing deep learning 2D / 3D image registration method is solved. Meanwhile, a pre-processing method based on Grangeat relationship and Radon domain transformation is proposed, and the robustness and generalization ability of the deep learning 2D / 3D image registration method are significantly improved.

Description

Technical Field

[0001] This invention relates to the field of image applications, and in particular to a deep learning-based 2D / 3D image registration method for pose matching of 2D X-ray images and 3D CT images. Background Technology

[0002] 2D / 3D image registration is an important image processing technique in image applications, applicable to various industrial fields requiring combined 2D / 3D image processing. For example, it can be used as a core component of robot-assisted surgical systems. In this field, 2D / 3D image registration is a technique for estimating the relative spatial pose between a known 3D structure and its corresponding 2D image. The 3D structure includes three-dimensional computed tomography (CT) images of skeletal anatomy or three-dimensional computer-aided drafting (CAD) models of surgical instruments, while the 2D image includes X-ray fluoroscopic images of the same 3D structure. Currently, deep learning-based 2D / 3D image registration methods directly predict the pose matrix corresponding to the X-ray image using a neural network model. The advantage of this method is its fast forward inference speed, significantly reducing registration time. However, existing deep learning-based methods require extensive retraining and adjustments for different types of images, resulting in high time and computational costs. The reason for this situation is that current deep learning registration methods are improvements based on iterative optimization of Digitally Reconstructed Radiograph (DRR) images. They utilize ray casting algorithms to project 2D DRR images from 3D CT data, directly using these projected DRR images as training images for the network and their corresponding pose parameters as labels to train a 2D / 3D image registration network model for the corresponding CT data. However, this model can only be used with that specific image data, requiring repeated operations for each different image dataset. This results in high training costs, long waiting times, and fails to meet practical application requirements. Furthermore, finding a suitable method to train both 3D and 2D images together is difficult, thus hindering the achievement of higher robustness and generalization capabilities in image processing. Summary of the Invention

[0003] Purpose of the invention: To address the above-mentioned shortcomings, this invention provides a 2D / 3D image registration method based on Radon domain transform, which trains a neural network to improve the robustness and generalization of deep learning-based 2D / 3D image registration methods.

[0004] Technical Solution: To solve the above problems, this invention discloses a 2D / 3D image registration method based on Radon domain transform, characterized by the following steps:

[0005] Step 1: Acquire the 3D CT volume data and 2D X-ray image of the target area;

[0006] Step 2: Preprocess the 3D CT volume data and 2D X-ray image data respectively; among which,

[0007] Preprocessing of 3D CT volume data involves transforming the image coordinate system of the 3D CT volume data to a 3D coordinate system and aligning the AP view projected from the positive Y-axis to the negative Y-axis with the direction of the 2D X-ray image, thus obtaining the processed CT position image data. TruePose ;

[0008] The preprocessing of a two-dimensional X-ray image involves: normalizing the pixel values ​​of the two-dimensional X-ray image to the range of [0,1], downsampling the two-dimensional X-ray image to a specified size, and then standardizing the pixel values ​​using a normal distribution with a given mean and standard deviation to obtain the preprocessed two-dimensional X-ray image I. xray ;

[0009] Step 3: Obtain several existing CT volume data of the same target region of the same type of object as the object in step (1) from the existing CT dataset as a training set, and preprocess the existing CT volume data, including:

[0010] Step 3.1: Obtain CT volume data - index group CT indexed ;

[0011] Step 3.2: Sample the true viewpoint pose parameters T from a normal distribution with given mean and standard deviation. gt Through this true viewpoint pose parameter T gt DRR image I is generated using a ray casting algorithm. gt Integrate all DRR images gt and the corresponding true viewpoint pose parameters T gt Obtain the DRR dataset I drr_dataset ;

[0012] Step 3.3: Analyze the DRR dataset I obtained in Step 3.2. drr_dataset DRR image I gt Preprocessing: For DRR image I gt By normalizing the maximum and minimum pixel values ​​of each image, the image resolution is adjusted to a specified size for the DRR image I. gtThe pixel values ​​are standardized using a normal distribution with a given mean and standard deviation to obtain the preprocessed DRR dataset I. drr_dataset_pre ;

[0013] Step 4: Perform 3D Radon domain transformation on the existing CT volume data from Step 3, and then calculate the 3D training dataset, i.e., the Radon dataset, using Grangerat relations. 3D Dataset Radon 3D Each data point is Radon 3D_inpud ;

[0014] For the DRR image I that has been preprocessed in step 3.3 gt The 2D training dataset, namely the Radon dataset, is obtained by performing a 2D Radon domain transformation and then calculating using Grangerat relations. 2D Dataset Radon 2D Includes Radon 2D_input and its corresponding pose parameter T gt ;

[0015] Step 5: Construct a deep neural network model for 2D / 3D image registration and design the loss function used during network training, using the Radon dataset. 2D Randomly selected data Radon 2D_input and from the Radon dataset 3D Find Rado 2D_input Corresponding Radon 3D_input Radon 2D_input Radon 3D_input As input to the model, the model will use Radon 2D_input Radon 3D_input The encoding is the corresponding viewpoint pose parameters, and the model output is the predicted viewpoint pose parameters T. pred ;

[0016] Step 6: Transfer the correct CT location image data from Step 2 to the CT scanner. TruePose The CT is obtained by performing a 3DRadon domain transformation and then calculating using the Grangerat relation. TruePose_Radom For the real X-ray image I in step 2 xray Perform a 2DRadon domain transformation, and then calculate I using the Grangerat relation. xray_Radon Load the trained model and use CT. TruePose_Radon I xray_Radon Input the trained model, and the model predicts I. xray The viewpoint pose parameters T coarse T coarseThe results show the 2D / 3D image registration based on deep learning predictions.

[0017] Furthermore, in step 3.1, the existing CT volume data is preprocessed by dividing it into regions with CT values ​​< a HU, regions with a HU < CT value < b HU, and regions with CT values ​​> b HU, where HU is the unit of CT value, and a and b are preset values. An index is created for each region corresponding to each voxel and saved, resulting in the preprocessed CT volume data - index group CT. indexed .

[0018] Furthermore, in step 3.2, during each forward projection, the CT volume data-index group CT created in step 3.1 is used. indexed For regions in CT where the CT value is greater than b HU, the voxel values ​​are uniformly multiplied by the bone enhancement coefficient μ, where μ is randomly sampled from a uniform distribution within a given range.

[0019] Furthermore, in step 4, the physical size of the existing CT volume data is calculated to obtain the diagonal physical length D. CTdiagonal Determine the size D of the sampling plane. CTdiagonal *D CTdiagonal According to the physical length D of the diagonal CTdiagonal The discrete sampling interval of the Radon integration plane is determined, and the discrete sampling angle is calculated. The sampling coordinates in physical space for each integration plane are calculated, and then the sampling coordinates are converted to CT index coordinates. Sampling coordinate points outside the range are excluded. Then, sampling and accumulation are performed on the CT data to obtain the integration result. All sampling points are integrated to obtain the 3DRadon domain transformation result. The result is then processed again according to the Grangerat relation and standardized using a normal distribution with a given mean and standard deviation to obtain the preprocessed 3DRadon dataset Radon. 3D ;

[0020] to I drr_dataset_pre Each DRR image in gt Perform the following processing: Calculate the DRR image I gt The true physical size, and calculate the DRR image I gt The actual diagonal length is used as the sampling length D DRRdiagonal According to the sampling length D DRRdiagonalThe discrete sampling interval of the Radon integration plane is determined, and the discrete sampling angle is calculated. The sampling coordinates in physical space for each integration plane are calculated, then converted to DRR image index coordinates. Sampling coordinates outside the range are excluded. The results are then sampled and accumulated in the DRR image to obtain the integration result. All sampling points are integrated to obtain the 2DRadon domain transform data. This data is then processed again according to the Grangerat relation and standardized using a normal distribution with given mean and standard deviation to obtain the preprocessed 2DRadon dataset (Radon). 2D .

[0021] Furthermore, in step 5, the deep neural network model comprises six sub-modules connected in a cascaded structure. The first module consists of a CBR layer and a 3*3*3 maxpool layer. The following four modules share the same minimum unit, differing only in the number of minimum units and the feature extraction depth. Each minimum unit consists of two CBR layers, a Conv layer, a BatchNorm layer, and residual connections. The last module is a fully connected layer composed of an AvgPool layer and a Linear layer, whose function is to regress the viewpoint pose parameter T from the high-dimensional features extracted from the first five sub-modules. pred .

[0022] Furthermore, the loss function used when constructing the network for training includes two functions.

[0023] The first loss function is the input Radon. 2D_input The corresponding label truth viewpoint pose parameter T gt and the viewpoint pose parameters T predicted by the model pred Position offset MSE loss

[0024]

[0025] Where dx gt For the current Radon 2D_input The corresponding offset of the CT pose relative to the center point X-axis, dx pred dy represents the X-axis offset predicted by the model. gt For the current Radon 2D_input The corresponding offset of the CT pose relative to the center point Y-axis, dy pred The Y-axis offset predicted by the model; dz gt For the current Radon 2D_input The corresponding offset of the CT pose relative to the center point Z-axis, dz pred This represents the Z-axis offset predicted by the model.

[0026] The second loss function is the input Radon.2D_input The corresponding label truth viewpoint pose parameter T gt and the viewpoint pose parameters T predicted by the model pred Angular offset MSE loss

[0027]

[0028] Among them yaw gt For the current Radon 2D_input The difference between the corresponding CT pose and the initial CT yaw angle, yaw pred The pitch is the difference in yaw angle predicted by the model. gt For the current Radon 2D_input The difference between the corresponding CT pose and the initial CT pitch angle, pitch pred The pitch angle difference predicted by the model; roll gt For the current Radon 2D_input The difference between the corresponding CT pose and the initial CT roll angle, roll pred This represents the difference in roll angle predicted by the model.

[0029] Furthermore, the model outputs the predicted viewpoint pose parameters T. pred It is a 1×6 matrix, the first three parameters are the spatial offset of CT relative to the center point, and the last three parameters are the rotational offset of CT from the initial pose to the current pose.

[0030] Furthermore, in step 5.3, the model is trained by backpropagation using the loss function constructed in step 5.2, and the model parameters are updated. The 2D / 3D image registration and prediction model is trained for a total of num_epoch epochs, with each epoch containing num_itr iterations. The model performs backpropagation to update parameters once in each iteration, and the model file is saved after training is completed.

[0031] Furthermore, the 2D / 3D image registration method based on Radon domain transform is used to predict the 2D / 3D image registration results through a deep neural network model. It is not used as a diagnostic criterion and is for non-diagnostic and non-therapeutic purposes.

[0032] Furthermore, the present invention also provides a 2D / 3D image registration system based on Radon domain transform, comprising the following steps:

[0033] The data acquisition module is used to acquire three-dimensional CT volume data and two-dimensional X-ray images of the target area of ​​the object;

[0034] Preprocessing module one is used to preprocess the 3D CT volume data and the 2D X-ray image data respectively; among them,

[0035] Preprocessing of 3D CT volume data involves transforming the image coordinate system of the 3D CT volume data to a 3D coordinate system and aligning the AP view projected from the positive Y-axis to the negative Y-axis with the direction of the 2D X-ray image, thus obtaining the processed CT position image data. TruePose ;

[0036] The preprocessing of a two-dimensional X-ray image involves: normalizing the pixel values ​​of the two-dimensional X-ray image to the range [0,1], downsampling the two-dimensional X-ray image to a specified size, and then standardizing the pixel values ​​using a normal distribution with a given mean and standard deviation to obtain the preprocessed two-dimensional X-ray image T. xray ;

[0037] Preprocessing module two is used to obtain several existing CT volume data of the same target region of the same type of object as the object in step 1 from the existing CT dataset as a training set, and to preprocess the existing CT volume data, including:

[0038] 3.1 Obtaining CT volume data - Index group CT indexed ;

[0039] 3.2 Obtaining the true viewpoint pose parameters T by sampling from a normal distribution with given mean and standard deviation. gt Through this true viewpoint pose parameter T gt DRR image I is generated using a ray casting algorithm. gt Integrate all DRR images gt and the corresponding true viewpoint pose parameters T gt Obtain the DRR dataset I drr_dataset ;

[0040] 3.3 For the acquired DRR dataset I drr_dataset DRR image I gt Preprocessing: For DRR image I gt By normalizing the maximum and minimum pixel values ​​of each image, the image resolution is adjusted to a specified size for the DRR image I. gt The pixel values ​​are standardized using a normal distribution with a given mean and standard deviation to obtain the preprocessed DRR dataset I. drr_dataset_pre ;

[0041] The domain transformation processing module performs 3DRadon domain transformation on existing CT volume data, and then calculates the 3D training dataset, i.e., the Radon dataset, using Grangerat relations. 3D Dataset Radon 3D Each data point is Radon 3D_input ;

[0042] For the preprocessed DRR image I gt The 2D training dataset, namely the Radon dataset, is obtained by performing a 2D Radon domain transformation and then calculating using Grangerat relations. 2D Dataset Radon 2D Includes Radon 2D_input and its corresponding pose parameter T gt ;

[0043] A deep neural network model for 2D / 3D image registration and a loss function module used during network training, derived from the Radon dataset. 2D Randomly selected data Radon 2D_input and from the Radon dataset 3D Find Radon 2D_input Corresponding Radon 3D_input Radon 2D_input Radon 3D_input As input to the model, the model will use Radon 2D_input Radon 3D_input The encoding is the corresponding viewpoint pose parameters, and the model output is the predicted viewpoint pose parameters T. pred ;

[0044] The model training and output module is used to train and output CT images of the correct CT location. TruePose The CT is obtained by performing a 3DRadon domain transformation and then calculating using the Grangerat relation. TruePose_Radon For real X-ray images I xray Perform a 2DRadon domain transformation, and then calculate I using the Grangerat relation. xTay_Radon Load the trained model and use CT. TruePose_Radon I xray_Radon Input the trained model, and the model predicts I. xray The viewpoint pose parameters T coarse T coarse The results show the 2D / 3D image registration based on deep learning predictions.

[0045] Beneficial effects: This invention addresses the problem of poor versatility in existing deep learning 2D / 3D image registration methods by proposing a model architecture that incorporates 3D structural information into the network input for network training. It implements a preprocessing method based on Grangerat relations and Radon domain transform, improving the robustness and generalization of deep learning-based 2D / 3D image registration methods. Attached Figure Description

[0046] Figure 1This is a schematic diagram of the overall process of the multi-stage 2D / 3D image registration method of the present invention.

[0047] Figure 2 (a) is a schematic diagram of the three-dimensional display effect of the acquired CT body data.

[0048] Figure 2 (b) is an unprocessed X-ray image.

[0049] Figure 2 (c) is a pre-processed X-ray image.

[0050] Figure 2 (d) is a preprocessed DRR image that has been registered with the X-ray image.

[0051] Figure 3 (a) is a real X-ray image.

[0052] Figure 3 (b) shows the incorrect initial CT pose and the AP view projected from the positive Y-axis direction to the negative Y-axis direction.

[0053] Figure 3 (c) shows the correct initial CT pose and the AP view projected from the positive Y-axis direction to the negative Y-axis direction.

[0054] Figure 4 (a) is the dataset Radon 3D A slice of data in the image.

[0055] Figure 4 (b) is the dataset Radon 2D An image from [the source].

[0056] Figure 5 This is a schematic diagram of the coordinate axes of a CT slice and the coordinate axes of a loading matrix.

[0057] Figure 6 This is a schematic diagram of the 2DRadon domain transformation process.

[0058] Figure 7 This is a schematic diagram of the 3DRadon domain transformation process.

[0059] Figure 8 This is a diagram illustrating the Grangate relationship.

[0060] Figure 9 This is a schematic diagram of the 2D / 3D image registration deep neural network architecture proposed in the technical solution of this invention. Detailed Implementation

[0061] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to embodiments of the human pelvic structure and accompanying drawings. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0062] The specific implementation process of the present invention is as follows: Figure 1 As shown, the steps are as follows:

[0063] Step 1: Obtain the 3D CT volume data and 2D X-ray image data of the target area. An example of 3D CT volume data is shown below. Figure 2 As shown in (a), an example of the effect of a single image of a two-dimensional X-ray is as follows. Figure 2 As shown in (b).

[0064] Step 2: Preprocess the CT data and X-ray image data respectively;

[0065] Step 2.1: Preprocess the 3D CT volume data, performing operations such as rotation and axis flipping on the CT. For example... Figure 5 As shown, the CT slice coordinate axes differ from the loading matrix coordinate axes. To ensure the accuracy of projection and Radon domain transformation, the slice coordinate axes must be aligned with the matrix coordinate axes. Simultaneously, it is also necessary to ensure that the AP view projected from the positive Y-axis to the negative Y-axis by the CT scan is aligned with the direction of the 2D X-ray image. An example of a single 2D X-ray image is shown below. Figure 3 As shown in (a), the erroneous CT location and the projected DRR image are as follows: Figure 3 (b) shows the processed CT scan. TruePose and the projected DRR image, such as Figure 3 As shown in (c).

[0066] Step 2.2: Preprocess the acquired two-dimensional X-ray image. The preprocessing is performed in the following order: normalize the image pixel values ​​to the range [0,1]; then perform trilinear interpolation downsampling on the X-ray image to a specified size (resolution 500×500 pixels); then standardize the pixel values ​​using a normal distribution with a given mean and standard deviation (mean = 0.3080, standard deviation = 0.1494) to obtain the preprocessed two-dimensional X-ray image I. xray I xray An image effect as Figure 2 As shown in (c).

[0067] Step 3: Preprocess the training set of 3D CT volume data by using Siddon's Ray-tracing Algorithm to process the CT data. indexedPerform forward projection to generate DRR images for deep learning training, and preprocess the DRR dataset;

[0068] Step 3.1: Preprocess the 3D CT volumetric data of the training set. To ensure the robustness of the network training results, the CT data needs to be augmented. Operation 1: Slice the CT data by depth and rotate each slice image; Operation 2: Slice the CT data by depth and flip each slice image along the X-axis. Combining these two operations yields 8 times the original CT data. The augmented CT volumetric data is then divided into three regions based on the CT value (unit: HU) of each voxel: air (CT value < -800 HU), soft tissue (-800 HU < CT value < 350 HU), and bone (CT value > 350 HU). An index is created for each voxel's corresponding region and saved, resulting in the preprocessed CT volumetric data - indexed CT set. indexed .

[0069] Step 3.2: Use the ray projection algorithm to analyze the CT scan. indexed A forward projection is performed to generate a DRR image, and the true viewpoint pose parameters T are obtained by sampling from a normal distribution with a given mean and standard deviation. gt T gt Composed of three rotational components and three translational components, the sampling code is as follows: torch.distributions.Normal(2,10) represents a normal distribution with a mean of 2 and a standard deviation of 10.

[0070] r1=torch.distributions.Normal(0,0.2).sample(),

[0071] r2=torch.distributions.Normal(0,0.1).sample(),

[0072] r3=torch.distributions.Normal(0,0.25).sample(),

[0073] t1=torch.distributions.Normal(10,70).sample(),

[0074] t2=torch.distributions.Normal(250,90).sample(),

[0075] t3=torch.distributions.Normal(5,50).sample().

[0076] Through this pose parameter T gt DRR image generated using ray casting algorithm gt Integrate all images gt and the corresponding pose parameters T gt Obtain the DRR dataset I drr_dataset Each augmented CT scan was performed 500 times, generating a total of 128*8*500 DRR images. Each DRR image has an I... gt With the corresponding viewpoint pose parameters T gt As labels, a total of 512,000 DRR image-viewpoint pose parameter labels are provided, as shown in the following formula.

[0077] T gt =[dx,dy,dz,yaw,pitch,roll]

[0078] In step 3.2, the CT volume data-index group CT created in step 3.1 is used during each forward projection. indexed The voxel values ​​of the skeletal region in the CT scan are uniformly multiplied by the bone enhancement coefficient μ. μ is randomly sampled from a uniform distribution in the range [1,10]. The DRR image is then projected from the CT scan after bone enhancement using a ray projection algorithm.

[0079] Step 3.3: Analyze the DRR dataset I obtained in Step 3.2. drr_dataset Image I gt Preprocessing is performed in the following order: The DRR images are normalized using the maximum and minimum pixel values ​​of each image; the image resolution is adjusted to a specified size (500×500 pixels); and the pixel values ​​of the DRR images are standardized using a normal distribution with a mean of 0.3080 and a standard deviation of 0.1494, resulting in the preprocessed DRR dataset I. drr_dataset_pre Dataset I drr_dataset_pre One of the images in the image is as follows Figure 2 As shown in (d).

[0080] Step 4: Perform 3DRadon domain transformation on the 3D CT data of the training set, and then calculate the 3D training dataset using Grangerat relations; perform 2DRadon domain transformation on the preprocessed DRR images from Step 3.3, and then calculate the 2D training dataset using Grangerat relations. The 2DRadon domain transformation process is as follows: Figure 6 As shown, the formula derivation is in Appendix 1; the 3DRadon domain transformation process is as follows: Figure 7 As shown, the formula derivation is in Appendix 2; the Grangerat relationship is as follows. Figure 8 As shown, the formula derivation is in Appendix 3;

[0081] Step 4.1: Calculate the true physical size of the CT (pixel dimension multiplied by the actual physical spacing between pixels), and calculate the true diagonal length D. CTdiagonal (length of the body diagonal) Determine the size D of the sampling plane CTdiagonal *D CTdiagonal According to D CTdiagonal The discrete sampling interval of the Radon integration plane is determined, and the integration plane is represented in polar coordinates (d, θ, ...). ), where d is the distance from the origin to the plane; θ is the angle between the plane normal and the Z-axis; The angle between the normal to the plane projected onto the XY plane and the X-axis. To ensure that the 3DRadon calculation results of all CT scans have the same shape, d is within the range Uniform sampling is performed within the range [0, π], where θ is uniformly sampled. Uniform sampling is performed within the range [0, 2π]. The physical sampling coordinates of each integral plane are calculated using the following formula, where SampleX and SampleY are discrete sampling points in the plane, within the range... Samples were taken at 1mm intervals.

[0082]

[0083] point_index = d*n + SampleX*e x +SampleY*e y

[0084] At this point, `point_index` represents the physical coordinates. Dividing it by the actual physical pixel spacing across the three dimensions of CT converts it to CT index coordinates. Sampling coordinates outside this range are excluded. CT is then sampled according to these index coordinates. The sampling code is `nn.functional.grid_sample(volume,volume_SamplingMatrix,mode='bilinear',padding_mode='zeros',align_corners=False)`. All sampling results are summed to obtain the integral result. All sampling points are integrated to obtain 3DRadon domain transformed data. The d-dimensional difference of the 3DRadon data is calculated using the Grangerat formula, and then standardized using a normal distribution with a mean of 0 and a standard deviation of 1. This yields the preprocessed 3DRadon dataset, Radon. 3D Dataset Radon 3D Each data point is Radon 3D_input Dataset Radon 3D One of the Radon 3D_input slices such as Figure 4 As shown in (a).

[0085] Step 4.2, for I drr_dataset_pre Each DRR image is processed as follows: Calculate the true physical size (pixel dimension multiplied by the actual physical spacing between pixels), and calculate the true diagonal length D. DRRdiagonal (face diagonal length, Determine the sampling length D DRRdiagonal According to D DRRdiagonal The discrete sampling interval of the Radon integration plane is determined. The two-dimensional Radon integration plane is a line integral, and the integration plane is represented in polar coordinates (d, θ), where d is the perpendicular distance from the origin to the integration line; θ is the angle between the normal to the integration line and the X-axis. To ensure that the 2DRadon calculation results of all DRR images have the same shape, d is within the range... Uniform sampling is performed within the integral plane, with θ uniformly sampled within the range [0, 2π]. The physical sampling coordinates of each integral plane are calculated using the following formula, where SampleT represents the discrete sampling points of the integral line within the range [0, 2π]. Samples were taken at 1mm intervals.

[0086] n = (cos(θ), sin(θ))

[0087] e = (sin(θ), -cos(θ))

[0088] point_index = d*n + SampleT*e

[0089] At this point, `point_index` represents the physical coordinates. Dividing it by the actual physical pixel spacing in both dimensions of the DRR image converts it to DRR image index coordinates. Sampling coordinates outside the range are excluded. The DRR image is then sampled according to the index coordinates. The sampling code is `nn.functional.grid_sample(drr_img,drr_img_SamplingMatrix,mode='bilinear',padding_mode='zeros',align_corners=False)`. All sampling results are summed to obtain the integral result. All sampling points are integrated to obtain 2DRadon domain transform data. The d-dimensional difference of the 2DRadon data is calculated using the Grangerat formula. Then, the cosine value of the angle between the plane formed by the current integration line and the source point and the line formed by the origin of the source and detector points, `cos(β)`, is calculated, multiplied by the square of its reciprocal, and finally standardized using a normal distribution with a mean of 0 and a standard deviation of 1. This yields the preprocessed 2DRadon dataset Radon. 2D Dataset Radon 2D Includes Radon 2D_inputand its corresponding pose parameter T gt Dataset Radon 2D One of the Radon 2D_input like Figure 4 As shown in (b).

[0090] Step 5: Construct a deep neural network model for 2D / 3D image registration, and design the loss function to be used during network training, using the created Radon dataset. 3D Radon 2D and viewpoint pose parameters T gt Labels are used to train the network;

[0091] Step 5.1: Construct a deep neural network model, ResNet3D-50, for 2D / 3D image registration. The model structure and the loss function used are as follows: Figure 9 As shown. From the Radon dataset 2D Randomly selected data Radon 2D_input and from the Radon dataset 3D Find Radon 2D_input Corresponding Radon 3D_input Radon 2D_input Radon 3D_input The model is used as input. It comprises six sub-modules connected in a cascaded structure. The first module consists of a CBR layer (Conv+BatchNorm+ReLU) and a 3*3*3 maxpool layer. The following four modules share the same minimum unit, differing only in the number of minimum units and the feature extraction depth. Their minimum units consist of two CBR layers, a Conv layer, a BatchNorm layer, and residual connections. The last module is a fully connected layer, composed of an AvgPool layer and a Linear layer, which regresses the viewpoint pose parameters T from the high-dimensional features extracted from the first five sub-modules. pred .

[0092] Step 5.2: Construct the two loss functions used for network training;

[0093] The first loss function is the input Radon. 2D_input The corresponding label truth viewpoint pose parameter T gt and the viewpoint pose parameters T predicted by the model pred Position offset MSE loss

[0094]

[0095] Where dx gt For the current Radon 2D_inputThe corresponding offset of the CT pose relative to the center point X-axis, dx pred dy represents the X-axis offset predicted by the model. gt For the current Radon 2D_input The corresponding offset of the CT pose relative to the center point Y-axis, dy pred The Y-axis offset predicted by the model; dz gt For the current Radon 2D_input The corresponding offset of the CT pose relative to the center point Z-axis, dz pred This is the Z-axis offset predicted by the model.

[0096] The second loss function is the input Radon. 2D_inpu The corresponding label truth viewpoint pose parameter T gt and the viewpoint pose parameters T predicted by the model pred Angular offset MSE loss

[0097]

[0098] Among them yaw gt For the current Radon 2D_input The difference between the corresponding CT pose and the initial CT yaw angle, yaw pred The pitch is the difference in yaw angle predicted by the model. gt For the current Radon 2D_input The difference between the corresponding CT pose and the initial CT pitch angle, pitch pred The pitch angle difference predicted by the model; roll gt For the current Radon 2D_input The difference between the corresponding CT pose and the initial CT roll angle, roll pred This represents the difference in roll angle predicted by the model.

[0099] Step 5.3: Use the created Radon dataset 3D Radon 2D and viewpoint pose parameters T gt Labels are used to train the network. The model's input is Radon. 3D_input Radon 2D_input The model encodes it into the corresponding viewpoint pose parameters, and the model output is the predicted viewpoint pose parameters T. pred , among which, T pred It is a 1×6 matrix, the first three parameters are the spatial offset of CT relative to the center point, and the last three parameters are the rotation offset of CT from the initial pose to the current pose.

[0100] In step 5.3, Figure 9The text shows the size of the output features for each submodule, where B represents the batch size used for network training. During model training, B is set to 20. The text also mentions Radon code as the model input. 3D_input Size is (20,4,49,91,181), Radon 2D_input The size (20, 4, 49, 181) indicates that the model uses 20 Radon in one iteration. 3D_input Radon 2D_input As input, the pose parameter matrix label for each input is shown in the example below:

[0101] T gt

[0102] = [-67.721634,426.78497,-34.029266,1.3629382,0.032742713,-1.5193983]

[0103] The feature size after the first submodule is (20, 64, 13, 23, 46), the feature size after the second submodule is (20, 256, 13, 23, 46), the feature size after the third submodule is (20, 512, 7, 12, 23), the feature size after the fourth submodule is (20, 1024, 4, 6, 12), the feature size after the fifth submodule is (20, 2048, 2, 3, 6), and the feature size after the sixth submodule is (20, 6), representing 20 rigid transformation pose parameter matrices T of size 1×6. pred T pred Here is an example.

[0104] T gt

[0105] = [-15.624178,238.30144,-71.173882,1.6189629,0.14988339,-1.4973321]

[0106] In step 5.3, the model is trained by backpropagation using the loss function constructed in step 5.2, and the model parameters are updated. The 2D / 3D image registration and prediction model is trained for a total of 20 epochs, with 25,600 iterations in each epoch. The Adam optimizer is used for optimization. The model performs backpropagation to update the parameters once in each iteration. After training is completed, the model is saved as a ".pth" file.

[0107] Step 6: CT scan TruePose The CT is obtained by performing a 3DRadon domain transformation and then calculating using the Grangerat relation. TruePose_RadonFor the preprocessed real X-ray image I xray Perform a 2DRadon domain transformation, and then calculate I using the Grangerat relation. xray_Radon Load the trained model and use CT. TruePose_Radon I xray_Radon Input the trained model, and the model predicts I. xray The viewpoint pose parameters T coarse T coarse T represents the 2D / 3D image registration results based on deep learning predictions. coarse Here is an example.

[0108] T coarse =

[0109] [-5.5778799,221.61359,-37.559978,1.5399648,0.034316484,-1.5442889]. Appendix 1: 2DRadon Domain Transformation Formula

[0110] The general form of the n-dimensional Radon domain transform is shown in the following vector integral:

[0111]

[0112] In the above formula Represents a point in n-dimensional space. It represents the unit normal vector of the (n-1)-dimensional hyperplane that is d equidistant from the origin.

[0113] In the two-dimensional Radon domain transformation, the integration hyperplane is a straight line.

[0114] make

[0115] n=(cos(θ),sin(θ))#(2)

[0116] The above formula can be transformed into:

[0117]

[0118] In the above formula, θ represents the angle between the unit normal vector and the X-axis in the two-dimensional Cartesian coordinate system, and (x,y) represents a point in the two-dimensional space in the two-dimensional Cartesian coordinate system.

[0119] Since the delta function is only valid when the expression is 0, the above surface integral can be transformed into a hyperplane linear integral, letting...

[0120] e=(sin(θ),-cos(θ))#(4)

[0121]

[0122] In the above formula, t represents a point on the integral line.

[0123] Appendix 2: Formula for 3DRadon Domain Transformation

[0124] The general form of the n-dimensional Radon domain transform is shown in the following vector integral:

[0125]

[0126] In the above formula Represents a point in n-dimensional space. It represents the unit normal vector of the (n-1)-dimensional hyperplane that is d equidistant from the origin.

[0127] In the three-dimensional Radon domain transformation, the integral hyperplane is a two-dimensional plane.

[0128] make

[0129]

[0130] The above formula can be transformed into:

[0131]

[0132] In the above formula, θ represents the angle between the unit normal vector and the Z-axis in the three-dimensional Cartesian coordinate system. The angle between the projection of the unit normal vector onto the XY plane in the three-dimensional Cartesian coordinate system and the X-axis is represented by (x,y,z), where (x,y,z) represents a point in the three-dimensional Cartesian coordinate system.

[0133] Since the delta function is only valid when the expression is 0, the above surface integral can be converted into a two-dimensional plane surface integral, let

[0134]

[0135] In the above formula, s x s y This represents a point on the integral plane.

[0136] Appendix 3: Grangate relation formula

[0137] Draw a straight line on the detector, and denote the integral of the cone beam projection data on the line as s(v).

[0138] Let the direction of the integration line on the planar detector be the u-axis direction. The direction perpendicular to it is the v-axis direction.

[0139] In the uv coordinate system, the cone beam data is denoted as g(u,v). If we represent the object f using a spherical coordinate system, the cone beam projection data g(u,v) can be expressed as:

[0140]

[0141] In the above formula, Let d represent the unit normal vector of the plane through which one ray of the cone beam passes, d represent the distance from the plane to the origin, D represent the distance of the line segment formed by the radiation source and the origin of the detector, θ represent the angle between the projection of the cone beam ray and the straight line formed by the radiation source and the origin of the detector onto the plane, and r represent a point on the cone beam ray.

[0142] Before calculating the data integral, all cone-beam projection data must be multiplied by a cosine function. Substituting into the above equation, the data integral on the u-axis is actually:

[0143]

[0144] Replacing variable u with variable θ, we can obtain formula 4 from formula 3:

[0145]

[0146] Establish the relationship between s(v) and the Radon domain transform of f. The Radon domain transform of f is:

[0147]

[0148] It also has the following relationships:

[0149]

[0150] In the above formula, β represents the angle between the straight line formed by the origin of the radiation source and the detector and the integration plane.

[0151] Thus we obtain:

[0152]

[0153] Finally, transform the partial derivative with respect to β into the partial derivative with respect to v:

[0154]

[0155] Ultimately, the relationship between the three-dimensional Radon domain transform and the two-dimensional Radon domain transform of volume data can be obtained:

[0156]

Claims

1. A 2D / 3D image registration method based on Radon domain transform, characterized in that, Includes the following steps: Step 1: Acquire the 3D CT volume data and 2D X-ray image of the target area; Step 2: Preprocess the 3D CT volume data and 2D X-ray image data respectively; among which, Preprocessing of 3D CT volume data involves transforming the image coordinate system of the 3D CT volume data to a 3D coordinate system and aligning the AP view projected from the positive Y-axis to the negative Y-axis with the orientation of the 2D X-ray image, thus obtaining the processed image data with the correct CT position. ; Preprocessing of a two-dimensional X-ray image involves: normalizing the pixel values ​​of the two-dimensional X-ray image to the range [0,1], downsampling the two-dimensional X-ray image to a specified size, and then standardizing the pixel values ​​using a normal distribution with a given mean and standard deviation to obtain the preprocessed two-dimensional X-ray image. ; Step 3: Obtain several existing CT volume data of the same target region of the same type of object as the object in step (1) from the existing CT dataset as a training set, and preprocess the existing CT volume data, including: Step 3.1: Obtain CT volume data - index group ; Step 3.2: Sample the true viewpoint pose parameters from a normal distribution with given mean and standard deviation. Through this true viewpoint pose parameter DRR images are generated using a ray casting algorithm. Integrate all DRR images and the corresponding true viewpoint pose parameters Obtain the DRR dataset ; Step 3.3: Process the DRR dataset obtained in Step 3.

2. DRR image in Preprocessing: DRR image By normalizing the maximum and minimum pixel values ​​of each image, the resolution of the image is adjusted to a specified size for the DRR image. The pixel values ​​are standardized using a normal distribution with a given mean and standard deviation to obtain the preprocessed DRR dataset. ; Step 4: Perform 3D DRadon domain transformation on the existing CT volume data from Step 3, and then calculate the 3D training dataset using Grangerat relations. Dataset Each data point is ; For the DRR image preprocessed in step 3.3 The 2D training dataset is obtained by performing a 2D DRadon domain transformation and then calculating using Grangerat relations. Dataset Include and its corresponding pose parameters ; Step 5: Construct a deep neural network model for 2D / 3D image registration and design the loss function used during network training, from the dataset. Randomly selected data and from the dataset turn up corresponding , , As input to the model, the model will , The encoding is the corresponding viewpoint pose parameters, and the model output is the predicted viewpoint pose parameters. ; Step 6: Transfer the correct CT location image data from Step 2. Perform a 3DRadon domain transformation, and then calculate using the Grangerat relation to obtain... For the real X-ray images in step 2 Perform a 2DRadon domain transformation, and then calculate using the Grangerat relation. Load the trained model and use , Input the trained model, and the model predicts... Viewpoint pose parameters , The results show the 2D / 3D image registration based on deep learning predictions.

2. The 2D / 3D image registration method based on Radon domain transform according to claim 1, characterized in that, In step 3.1, the existing CT volume data is preprocessed by dividing it into regions with CT values ​​< a HU, a HU < CT value < b HU, and CT values ​​> b HU, where HU is the unit of CT value, and a and b are preset values. An index is created for each region corresponding to each voxel and saved, resulting in the preprocessed CT volume data-index group. .

3. The 2D / 3D image registration method based on Radon domain transform according to claim 2, characterized in that, In step 3.2, the CT volume data-index group created in step 3.1 is used during each forward projection. For regions in CT scans with CT values ​​greater than bHU, the voxel values ​​are uniformly multiplied by the bone enhancement factor. , Random sampling is performed from a uniform distribution within a given range.

4. The 2D / 3D image registration method based on Radon domain transform according to claim 3, characterized in that, In step 4, the physical size of the existing CT volume data is calculated to obtain the diagonal physical length. Determine the size of the sampling plane Based on the physical length of the diagonal The discrete sampling interval of the Radon integration plane is determined, and the discrete sampling angle is calculated. The sampling coordinates in physical space for each integration plane are calculated, and then the sampling coordinates are converted to CT index coordinates. Sampling coordinate points outside the range are excluded. Then, sampling and accumulation are performed on the CT data to obtain the integration result. All sampling points are integrated to obtain the 3DRadon domain transformation result. The result is then processed again according to the Grangerat relation and standardized using a normal distribution with a given mean and standard deviation to obtain the preprocessed 3DRadon dataset. ; right Each DRR image Perform the following processing: Calculate the DRR image The true physical size, and calculate the DRR image. The actual diagonal length is used as the sampling length. According to the sampling length The discrete sampling interval of the Radon integration plane is determined, and the discrete sampling angle is calculated. The sampling coordinates in physical space for each integration plane are calculated, then converted to DRR image index coordinates. Sampling coordinates outside the range are excluded. The samples are then accumulated in the DRR image to obtain the integration result. All sampling points are integrated to obtain 2DRadon domain transform data. This data is then processed again according to the Grangerat relation and standardized using a normal distribution with given mean and standard deviation to obtain the preprocessed 2DRadon dataset. .

5. The 2D / 3D image registration method based on Radon domain transform according to claim 1, characterized in that, In step 5, the deep neural network model contains six sub-modules, which are connected in a series structure; the first module consists of a CBR layer and a... The first module consists of a maxpool layer; the next four modules have the same minimum unit, differing only in the number of minimum units and the feature extraction depth. Each minimum unit consists of two CBR layers, a Conv layer, a BatchNorm layer, and residual connections; the last module is a fully connected layer, consisting of an AvgPool layer and a Linear layer, whose function is to regress the viewpoint pose parameters from the high-dimensional features extracted from the first five sub-modules. .

6. The 2D / 3D image registration method based on Radon domain transform according to claim 5, characterized in that, The loss functions used when building the network for training include two types. The first loss function is the input. Corresponding label truth viewpoint pose parameters and the viewpoint pose parameters predicted by the model Position offset MSE loss , ， in For the present The corresponding offset of the CT pose relative to the X-axis of the center point. The X-axis offset predicted by the model; For the present The corresponding offset of the CT pose relative to the center point on the Y-axis. This represents the Y-axis offset predicted by the model. For the present The corresponding offset of the CT pose relative to the Z-axis of the center point. This represents the Z-axis offset predicted by the model. The second loss function is the input Corresponding label truth viewpoint pose parameters and the viewpoint pose parameters predicted by the model Angular offset MSE loss , ， in For the present The difference between the corresponding CT pose and the initial CT yaw angle. The difference in yaw angle predicted by the model; For the present The difference between the corresponding CT pose and the initial CT pitch angle. The pitch angle difference predicted by the model; For the present The difference between the corresponding CT pose and the initial CT roll angle. This represents the difference in roll angle predicted by the model.

7. The 2D / 3D image registration method based on Radon domain transform according to claim 6, characterized in that, Predicted viewpoint pose parameters output by the model It is a 1×6 matrix, the first three parameters are the spatial offset of CT relative to the center point, and the last three parameters are the rotational offset of CT from the initial pose to the current pose.

8. The 2D / 3D image registration method based on Radon domain transform according to claim 5, characterized in that, In step 5, the model is trained by backpropagation using the constructed loss function, and the model parameters are updated. The 2D / 3D image registration and prediction model is trained for a total of num_epoch epochs, with each epoch containing num_itr iterations. The model performs backpropagation to update the parameters once during each iteration. After training is complete, the model file is saved.

9. The 2D / 3D image registration method based on Radon domain transform according to any one of claims 1 to 8, characterized in that, The 2D / 3D image registration method based on Radon domain transform is used to predict 2D / 3D image registration results through a deep neural network model. It is not used as a diagnostic standard and is for non-diagnostic and non-therapeutic purposes.

10. A 2D / 3D image registration system based on Radon domain transform, characterized in that, Includes the following steps: The data acquisition module is used to acquire three-dimensional CT volume data and two-dimensional X-ray images of the target area of ​​the object; Preprocessing module one is used to preprocess the 3D CT volume data and the 2D X-ray image data respectively; among them, Preprocessing of 3D CT volume data involves transforming the image coordinate system of the 3D CT volume data to a 3D coordinate system and aligning the AP view projected from the positive Y-axis to the negative Y-axis with the orientation of the 2D X-ray image, thus obtaining the processed image data with the correct CT position. ; Preprocessing of a two-dimensional X-ray image involves: normalizing the pixel values ​​of the two-dimensional X-ray image to the range [0,1], downsampling the two-dimensional X-ray image to a specified size, and then standardizing the pixel values ​​using a normal distribution with a given mean and standard deviation to obtain the preprocessed two-dimensional X-ray image. ; Preprocessing module two is used to obtain several existing CT volume data of the same target region of the same type of object as the object in step 1 from the existing CT dataset as a training set, and to preprocess the existing CT volume data, including: 3.1 Obtaining CT volume data - index group ; 3.2 Obtaining True Viewpoint Pose Parameters by Sampling from a Normal Distribution with Given Mean and Standard Deviation Through this true viewpoint pose parameter DRR images are generated using a ray casting algorithm. Integrate all DRR images and the corresponding true viewpoint pose parameters Obtain the DRR dataset ; 3.3 On the acquired DRR dataset DRR image in Preprocessing: DRR image By normalizing the maximum and minimum pixel values ​​of each image, the resolution of the image is adjusted to a specified size for the DRR image. The pixel values ​​are standardized using a normal distribution with a given mean and standard deviation to obtain the preprocessed DRR dataset. ; The domain transformation processing module performs 3DRadon domain transformation on existing CT volume data, and then calculates the 3D training dataset using Grangerat relations. Dataset Each data point is ; For the preprocessed DRR image The 2D training dataset is obtained by performing a 2D DRadon domain transformation and then calculating using Grangerat relations. Dataset Include and its corresponding pose parameters ; A deep neural network model for 2D / 3D image registration and a loss function module used during network training to extract data from the dataset. Randomly selected data and from the dataset turn up corresponding , , As input to the model, the model will , The encoding is the corresponding viewpoint pose parameters, and the model output is the predicted viewpoint pose parameters. ; The model training and output module is used to train and output image data at the correct CT location. Perform a 3DRadon domain transformation, and then calculate using the Grangerat relation to obtain... For real X-ray images Perform a 2DRadon domain transformation, and then calculate using the Grangerat relation. Load the trained model and use , Input the trained model, and the model predicts... Viewpoint pose parameters , The results show the 2D / 3D image registration based on deep learning predictions.

Images