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

Abstract

The invention discloses a 2D / 3D image registration method and system based on Radon domain transformation. The registration method comprises two key components: a pretreatment technology based on Radon domain transformation and neural network model prediction. According to the method, 3D structure information is added into network input, and the network is trained, so that the problem of poor universality of an existing deep learning 2D / 3D image registration method is solved. Meanwhile, a preprocessing method based on the Grangeat relation and Radon domain transformation is provided, and the robustness and generalization ability of the deep learning 2D / 3D image registration method are remarkably improved.

Description

Technical Field

[0001] The present invention relates to the field of image application, and in particular to a deep learning 2D / 3D image registration method for posture matching of two-dimensional X-ray images and three-dimensional CT images. Background Art

[0002] 2D / 3D image registration technology is an important image processing technology in the field of image application, and can be applied to multiple industrial fields that require 2D / 3D image combined processing. For example, it can be applied to the core components of robot-assisted surgery systems. In this field, 2D / 3D image registration is a technology that estimates the relative spatial pose relationship between a known 3D structure and its corresponding 2D image. The 3D structure includes a 3D computed tomography (CT) image of a skeletal anatomical structure or a 3D computer-aided drafting (CAD) model of a surgical instrument, and the 2D image includes an X-ray fluoroscopic image of the same 3D structure. At present, the 2D / 3D image registration method based on deep learning directly predicts the corresponding pose matrix of the X-ray image through a neural network model. The advantage of this method is that the forward reasoning speed of the model is fast, which can greatly reduce the registration time. However, the existing deep learning-based methods require a lot of retraining and adjustment for different types of images, resulting in high time and computational costs. The reason for this situation is that the current deep learning registration method is improved on the basis of iterative optimization of digitally reconstructed radiograph (DRR). The ray casting algorithm is used to project the three-dimensional CT data to obtain a two-dimensional DRR image. The projected DRR image is directly used as the image for network training, and the corresponding posture parameters are used as the label for network training to train the 2D / 3D image registration network model corresponding to the CT data. However, this model can only be used for this image data, which will cause each different image data to repeat the above operations, resulting in huge training costs, long waiting time, and not meeting actual usage requirements. At the same time, it is difficult to find a suitable method for putting 3D images and 2D images into model training together, so it is difficult to meet the needs of higher robustness and generalization of image processing. Summary of the invention

[0003] Purpose of the invention: In view of the above shortcomings, the present 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 the 2D / 3D image registration method based on deep learning.

[0004] Technical solution: To solve the above problems, the present invention discloses a 2D / 3D image registration method based on Radon domain transformation, which is characterized by comprising the following steps:

[0005] Step 1, obtaining three-dimensional CT volume data and two-dimensional X-ray images of the target area of ​​the object;

[0006] Step 2: Preprocess the three-dimensional CT volume data and the two-dimensional X-ray image data respectively; wherein:

[0007] The preprocessing of the 3D CT volume data is as follows: the image coordinate system of the 3D CT volume data is converted to the 3D coordinate system and the AP view projected from the positive direction of the Y axis to the negative direction of the Y axis is consistent with the direction of the 2D X-ray image, so as to obtain the processed correct CT position image data CT TruePose ;

[0008] The preprocessing of the two-dimensional X-ray image is as follows: 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 ​​through a normal distribution with a given mean and standard deviation to obtain the preprocessed two-dimensional X-ray image I xray ;

[0009] Step 3: obtaining a number of existing CT volume data of the same target area of ​​the same type of object as the object in step (1) in the existing CT data set as a training set, and preprocessing 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 parameter T from a normal distribution with a given mean and standard deviation gt , through the true viewpoint pose parameter T gt , using the ray casting algorithm to generate DRR image I gt , integrating all DRR images I gt And the corresponding true viewpoint pose parameter T gt Get DRR dataset I drr_dataset ;

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

[0013] Step 4: Perform 3DRadon domain transformation on the existing CT volume data in step 3, and then calculate the 3D training data set by Grangeat relation, i.e., the data set Radon 3D ; Dataset Radon 3D Each data is Radon 3D_inpud ;

[0014] For the DRR image I preprocessed in step 3.3 gt Perform 2DRadon domain transformation, and then calculate the 2D training data set by Grangeat relation, that is, the data set Radon 2D ; Dataset Radon 2D Includes Radon 2D_input And its corresponding pose parameter T gt ;

[0015] Step 5: Build a deep neural network model for 2D / 3D image registration and the loss function used in network training. 2D Randomly select data from Radon 2D_input , and from the Radon dataset 3D Find Rado 2D_input Corresponding Radon 3D_input , Radon 2D_input , Radon 3D_input As model input, the model takes Radon 2D_input , Radon 3D_input Encoded as the corresponding viewpoint pose parameters, the model output is the predicted viewpoint pose parameter T pred ;

[0016] Step 6: The correct CT position image data CT in step 2 TruePose Perform 3DRadon domain transformation and then calculate CT using the Grangeat relation TruePose_Radom , for the real X-ray image I in step 2 xray Perform 2DRadon domain transformation and then calculate I using the Grangeat relation xray_Radon , load the trained model and use CT TruePose_Radon ,I xray_Radon Input the trained model, the model predicts I xray The viewpoint pose parameter T coarse , T coarse2D / 3D image registration results based on deep learning prediction.

[0017] Further, in step 3.1, the existing CT volume data is preprocessed, that is, the existing CT volume data is preprocessed and divided into a region with CT value < a HU, a region with a HU < CT value < b HU, and a region with CT value > b HU according to the CT value of each voxel, where HU is the unit of CT value, a and b are both preset values, and an index of the region corresponding to each voxel is made and saved to obtain the preprocessed CT volume data - index group CT indexed .

[0018] Further, in step 3.2, at each forward projection, the CT volume data-index group CT generated in step 3.1 is used. indexed The voxel values ​​in the region with CT values ​​> b HU in CT are uniformly multiplied by the bone enhancement coefficient μ, which 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 physical length D of the diagonal line. CTdiagonal , determine the sampling plane size D CTdiagonal *D CTdiagonal ; According to the physical length of the diagonal D CTdiagonal Determine the discrete sampling spacing of the Radon integration plane, calculate the discrete sampling angle, calculate the sampling coordinates in the physical space of each integration plane, and then convert the sampling coordinates to CT index coordinates, exclude the sampling coordinate points that are out of range, and then accumulate the samples in the CT data to get the integration results, integrate all the sampling points to get the 3DRadon domain transformation results, process them again according to the Grangeat relationship, and standardize them through the normal distribution with a given mean and standard deviation to get the preprocessed 3DRadon data set Radon 3D ;

[0020] to I drr_dataset_pre Each DRR image I gt The following processing is performed: Calculate the DRR image I gt The real physical size of the image and calculate the DRR image I gt The actual diagonal length is taken as the sampling length D DRRdiagonal ; According to the sampling length D DRRdiagonalDetermine the discrete sampling spacing of the Radon integral plane and calculate the discrete sampling angle; calculate the sampling coordinates in the physical space of each integral plane, and then convert the sampling coordinates into DRR image index coordinates, exclude the sampling coordinate points that are out of range, and then accumulate the samples in the DRR image to obtain the integral result, integrate all the sampling points to obtain the 2DRadon domain transformation data, process it again according to the Grangeat relationship, and standardize it through a normal distribution with a given mean and standard deviation to obtain the preprocessed 2DRadon data set Radon 2D .

[0021] Furthermore, in step 5, the deep neural network model contains six submodules, which are connected in series. The first module consists of a CBR layer and a 3*3*3 maxpool layer. The next four modules have the same minimum unit, which differs only in the number of minimum units and feature extraction depth of each module. The minimum unit consists of two CBR layers plus a Conv layer and a BatchNorm layer and residual connection. The last module is a fully connected layer, which consists of an AvgPool layer and a Linear layer. Its function is to regress the viewpoint pose parameter T from the high-dimensional features extracted from the first five submodules. pred .

[0022] Furthermore, there are two loss functions used in network training:

[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 shift MSE loss

[0024]

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

[0026] The second loss function is the input Radon2D_input The corresponding label truth viewpoint pose parameter T gt And the viewpoint pose parameters T predicted by the model pred Angle 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 yaw angle difference predicted by the model; pitch gt For the current Radon 2D_input The difference between the corresponding CT posture and the initial CT pitch angle, pitch pred is the pitch angle difference predicted by the model; roll gt For the current Radon 2D_input The difference between the corresponding CT posture and the initial CT roll angle, roll pred is the roll angle difference predicted by the model.

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

[0030] Furthermore, in the step 5.3, the model is back-propagated and trained using the loss function constructed in step 5.2 to update the model parameters. The 2D / 3D image registration prediction model is trained for a total of num_epoch epochs, each epoch contains num_itr iterations, and the model is back-propagated and updated for each iteration. The model file is saved after the 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] In addition, the present invention also provides a 2D / 3D image registration system based on Radon domain transformation, comprising the following steps:

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

[0034] The preprocessing module 1 is used to preprocess the three-dimensional CT volume data and the two-dimensional X-ray image data respectively; wherein,

[0035] The preprocessing of the 3D CT volume data is as follows: the image coordinate system of the 3D CT volume data is converted to the 3D coordinate system and the AP view projected from the positive direction of the Y axis to the negative direction of the Y axis is consistent with the direction of the 2D X-ray image, so as to obtain the processed correct CT position image data CT TruePose ;

[0036] The preprocessing of the two-dimensional X-ray image is as follows: 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 ​​through a normal distribution with a given mean and standard deviation to obtain the preprocessed two-dimensional X-ray image T xray ;

[0037] The second preprocessing module is used to obtain a number of existing CT volume data of the same target area of ​​the same type of object as the object in step 1 in the existing CT data set as a training set, and preprocess the existing CT volume data, including:

[0038] 3.1 Obtaining CT Volume Data - Index Group CT indexed ;

[0039] 3.2 Sampling the true viewpoint pose parameter T from a normal distribution with a given mean and standard deviation gt , through the true viewpoint pose parameter T gt , using the ray casting algorithm to generate DRR image I gt , integrating all DRR images I gt And the corresponding true viewpoint pose parameter T gt Get DRR dataset I drr_dataset ;

[0040] 3.3 Acquired DRR Dataset I drr_dataset DRR image I gt Preprocessing: DRR image I gt By normalizing the maximum and minimum pixel values ​​of each image, the resolution of the image is adjusted to the specified size, and the DRR image I gt The pixel values ​​are standardized by 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 is used to perform 3DRadon domain transformation on the existing CT volume data, and then calculate the 3D training data set by Grangeat relationship, that is, the data set Radon 3D ; Dataset Radon 3D Each data is Radon 3D_input ;

[0042] For the pre-processed DRR image I gt Perform 2DRadon domain transformation, and then calculate the 2D training data set by Grangeat relation, that is, the data set Radon 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 in network training to extract the Radon dataset 2D Randomly select data from Radon 2D_input , and from the Radon dataset 3D Find Radon 2D_input Corresponding Radon 3D_input , Radon 2D_input , Radon 3D_input As model input, the model takes Radon 2D_input , Radon 3D_input Encoded as the corresponding viewpoint pose parameters, the model output is the predicted viewpoint pose parameter T pred ;

[0044] Model training and output module, used to convert the correct CT position image data CT TruePose Perform 3DRadon domain transformation and then calculate CT using the Grangeat relation TruePose_Radon , for the real X-ray image I xray Perform 2DRadon domain transformation and then calculate I using the Grangeat relation xTay_Radon , load the trained model and use CT TruePose_Radon ,I xray_Radon Input the trained model, the model predicts I xray The viewpoint pose parameter T coarse , T coarse 2D / 3D image registration results based on deep learning prediction.

[0045] Beneficial effects: In view of the poor versatility of existing deep learning 2D / 3D image registration methods, the present invention proposes a model architecture for training the network by adding 3D structural information to the network input. A pre-processing method based on Grangeat relationship and Radon domain transformation is implemented, which improves the robustness and generalization of deep learning-based 2D / 3D image registration methods. BRIEF DESCRIPTION OF THE DRAWINGS

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

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

[0048] Figure 2 (b) is the X-ray image without preprocessing.

[0049] Figure 2 (c) is the preprocessed X-ray image.

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

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

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

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

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

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

[0056] Figure 5 It is a schematic diagram of the CT slice coordinate axis direction and the loading matrix coordinate axis direction.

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

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

[0059] Figure 8 This is a diagram of the Grangeat relationship.

[0060] Fig. 9 It is a schematic diagram of the 2D / 3D image registration deep neural network architecture proposed in the technical solution of the present invention. DETAILED DESCRIPTION

[0061] In order to make the purpose, technical solution and advantages of the present invention clearer, the present invention is further described in detail below in conjunction with the human pelvic structure object embodiment and the accompanying drawings. It should be understood that the specific embodiments described herein are only for explaining the present invention and are not intended to limit the present invention.

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

[0063] Step 1: Obtain three-dimensional CT volume data and two-dimensional X-ray image data of the target area of ​​the object. An example of three-dimensional CT volume data is as follows: Figure 2 (a) shows an example of a two-dimensional X-ray image. Figure 2 (b) as shown.

[0064] Step 2, preprocessing the CT data and the X-ray image data respectively;

[0065] Step 2.1: Preprocess the 3D CT volume data by rotating the CT and performing axis flipping operations. Figure 5 As shown in the figure, the CT slice coordinate axis direction is different from the loading matrix coordinate axis direction. To ensure the accuracy of projection and Radon domain transformation, the slice coordinate axis and the matrix coordinate axis direction need to be consistent. At the same time, it is also necessary to ensure that the AP view projected from the positive direction of the Y axis to the negative direction of the Y axis is consistent with the direction of the two-dimensional X-ray image. An example of a two-dimensional X-ray image effect is shown below. Figure 3 As shown in (a), the wrong CT position and the projected DRR image are as follows Figure 3 (b) shows that the processed CT TruePose And the projected DRR image is as follows Figure 3 (c) as shown.

[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 of [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 like Figure 2 (c) as shown.

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

[0068] Step 3.1. Preprocess the three-dimensional CT volume data of the training set. In order to ensure that the results learned by the network training have high robustness, the CT data needs to be expanded. Operation 1: Slice the CT by depth and rotate each slice image; Operation 2: Slice the CT by depth and flip the data of each slice image along the X-axis. The two combinations can obtain 8 times the original CT data. The expanded CT volume data is divided into three regions according to the CT value of each voxel (unit: HU): air (CT value <-800HU), soft tissue (-800HU <CT value <350HU) and bone (CT value >350HU). Make an index for the corresponding area of ​​each voxel and save it to obtain the preprocessed CT volume data - index group CT indexed .

[0069] Step 3.2: Use the ray casting algorithm to CT indexed Perform forward projection to generate the DRR image, and sample the true viewpoint pose parameter T from a normal distribution with a given mean and standard deviation gt , T gt It consists of three rotation components and three translation 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 the pose parameter T gt Generating DRR images using ray casting algorithm I gt , integrating all images I gt And the corresponding pose parameter T gt Get DRR dataset I drr_dataset The expanded CT images are sampled 500 times each, generating a total of 128*8*500 DRR images. gt With the corresponding viewpoint pose parameter T gt As labels, there are a total of 512,000 sets of DRR image-viewpoint pose parameter labels, and the labels are 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 generated in step 3.1 is used for each forward projection. indexed The voxel values ​​of the bone region in the CT are uniformly multiplied by the bone enhancement coefficient μ, which is randomly sampled from a uniform distribution in the range of [1,10]. The DRR image is projected from the CT after bone enhancement by a ray casting algorithm.

[0079] Step 3.3: 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 by the maximum and minimum pixel values ​​of each image, the image resolution is adjusted to the specified size (resolution is 500×500 pixels), and the pixel values ​​of the DRR images are standardized by a normal distribution with a mean of 0.3080 and a standard deviation of 0.1494 to obtain the preprocessed DRR dataset I drr_dataset_pre , Dataset I drr_dataset_pre An image in Figure 2 (d) as shown.

[0080] Step 4: Perform 3DRadon domain transformation on the three-dimensional CT data of the training set, and then calculate the 3D training data set by Grangeat relation; perform 2DRadon domain transformation on the DRR image preprocessed in step 3.3, and then calculate the 2D training data set by Grangeat relation. Figure 6 As shown, the formula derivation is shown in Appendix 1; the 3DRadon domain transformation process is as follows Figure 7 As shown, the formula derivation is shown in Appendix 2; Grangeat relationship is as follows Figure 8 As shown, the formula derivation is shown in Appendix 3;

[0081] Step 4.1. Calculate the real physical size of the CT (pixel dimension multiplied by the actual physical spacing of the pixels) and calculate the real diagonal length D CTdiagonal (body diagonal length, Determine the sampling plane size D CTdiagonal *D CTdiagonal According to D CTdiagonal Determine the discrete sampling spacing of the Radon integration plane, which is represented by 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; is the angle between the projection of the plane normal to the XY plane and the X axis. To ensure that the 3DRadon calculation results of all CTs have the same shape, d is in the range uniformly sampled in the range [0, π], θ is uniformly sampled in the range [0, π], Sample uniformly in the range [0, 2π]. Calculate the physical sampling coordinates of each integration plane using the following formula, where SampleX and SampleY are discrete sampling points in the plane, in the range The samples were collected at 1mm intervals.

[0082]

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

[0084] At this time, point_index is the coordinate in the physical space. Divide it by the actual physical spacing of the pixels in the three dimensions of CT to convert it into the CT index coordinate. Exclude the sampling coordinate points that are out of range, and sample CT according to the index coordinate. The sampling code is nn.functional.grid_sample(volume, volume_SamplingMatrix, mode = 'bilinear', padding_mode = 'zeros', align_corners = False). Add all sampling results to get the integral result. Integrate all sampling points to get the 3DRadon domain transformation data, perform differential calculation on the d dimension of the 3DRadon data according to the Grangeat formula, and then standardize it through a normal distribution with a mean of 0 and a standard deviation of 1. Get the preprocessed 3DRadon dataset Radon 3D , Dataset Radon 3D Each data is Radon 3D_input , Dataset Radon 3D A Radon 3D_input The slices are as follows Figure 4 (a) shown.

[0085] Step 4.2: I drr_dataset_pre Each DRR image in is processed as follows: Calculate the real physical size (pixel dimension size multiplied by the actual physical spacing of pixels) and calculate the real diagonal length D DRRdiagonal (Diagonal length of the face, Determine the sampling length D DRRdiagonal According to D DRRdiagonal Determine the discrete sampling spacing of the Radon integration plane. The two-dimensional Radon integration plane is a line integral. The integration plane is represented by polar coordinates (d, θ), where d is the vertical distance from the origin to the integration line; θ is the angle between the normal of the integration line and the X-axis. To ensure that the 2DRadon calculation results of all DRR images have the same shape, d is in the range The physical sampling coordinates of each integration plane are calculated as follows, where SampleT is the discrete sampling point of the integration line, in the range The samples were collected at 1mm intervals.

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

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

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

[0089] At this time, point_index is the coordinate in the physical space. Divide it by the actual physical spacing of the pixels in the two dimensions of the DRR image to convert it into the index coordinate of the DRR image. Exclude the sampling coordinate points that are out of range, sample the DRR image according to the index coordinate. The sampling code is nn.functional.grid_sample(drr_img,drr_img_SamplingMatrix,mode='bilinear',padding_mode='zeros',align_corners=False), and add all the sampling results to get the integral result. Integrate all sampling points to get the 2DRadon domain transformation data, perform differential calculation on the d dimension of the 2DRadon data according to the Grangeat formula, and then calculate the cosine value cos(β) of the angle between the current integral line and the plane formed by the source point and the line formed by the origin of the source point detector, and then multiply it by the square of its reciprocal, and finally standardize it through a normal distribution with a mean of 0 and a standard deviation of 1. Get the preprocessed 2DRadon dataset Radon 2D , Dataset Radon 2D Includes Radon 2D_inputAnd its corresponding pose parameter T gt , Dataset Radon 2D A Radon 2D_input like Figure 4 (b) as shown.

[0090] Step 5: Build a deep neural network model for 2D / 3D image registration, build the loss function used in network training, and use the Radon dataset 3D , Radon 2D and viewpoint pose parameter T gt Labels are used to train the network;

[0091] Step 5.1: Build a deep neural network model resnet3D-50 for 2D / 3D image registration. The model structure and loss function used are as follows: Fig. 9 As shown. From the Radon dataset 2D Randomly select data from Radon 2D_input , and from the Radon dataset 3D Find Radon 2D_input Corresponding Radon 3D_input , Radon 2D_input , Radon 3D_input As model input. The model contains six sub-modules, which are connected in series. The first module consists of a CBR layer (Conv+BatchNorm+Relu) and a 3*3*3 maxpool layer. The next four modules have the same minimum unit, the only difference is the number of minimum units in each module and the feature extraction depth. The minimum unit consists of two CBR layers plus a Conv layer and a BatchNorm layer and a residual connection. The last module is a fully connected layer, which consists of an AvgPool layer and a Linear layer. Its function is to regress the viewpoint pose parameter T from the high-dimensional features extracted by the first five sub-modules. pred .

[0092] Step 5.2: Construct two loss functions 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 shift MSE loss

[0094]

[0095] where dx gt For the current Radon 2D_inputThe corresponding CT pose is offset from the X axis of the center point, dx pred The X-axis offset predicted by the model; dy gt For the current Radon 2D_input The corresponding CT pose is offset from the Y axis of the center point, dy pred The Y-axis offset predicted by the model; dz gt For the current Radon 2D_input The corresponding CT pose offset relative to the Z axis of the center point, dz pred 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 Angle 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 yaw angle difference predicted by the model; pitch gt For the current Radon 2D_input The difference between the corresponding CT posture and the initial CT pitch angle, pitch pred is the pitch angle difference predicted by the model; roll gt For the current Radon 2D_input The difference between the corresponding CT posture and the initial CT roll angle, roll pred is the roll angle difference predicted by the model.

[0099] Step 5.3: Use the Radon dataset 3D , Radon 2D and viewpoint pose parameter T gt The input of the model 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 parameter T pred , where T pred It is a 1×6 matrix. The first three parameters are the spatial offset of CT compared to the center point, and the last three parameters are the rotation offset of CT from the initial posture to the current posture.

[0100] In step 5.3, Fig. 9The output feature size of each submodule is shown in Figure 2. B represents the batch size used for network training. B is set to 20 during model training. 3D_input Size is (20,4,49,91,181), Radon 2D_input The size is (20,4,49,181), which means that the model uses 20 Radons in one iteration. 3D_input , Radon 2D_input As input, an example of the pose parameter matrix label corresponding to each input is:

[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) indicating 20 rigid transformation pose parameter matrices T of size 1×6 pred , T pred An example of is as follows.

[0104] T gt

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

[0106] In the step 5.3, the model is back-propagated and trained using the loss function constructed in step 5.2 to update the model parameters. The 2D / 3D image registration prediction model is trained for a total of 20 epochs, each epoch contains 25,600 iterations, and is optimized using the Adam optimizer. The model is back-propagated and updates the parameters once each iteration. After training, the model ".pth" file is saved.

[0107] Step 6: CT TruePose Perform 3DRadon domain transformation and then calculate CT using the Grangeat relation TruePose_Radon, for the real X-ray image I after preprocessing xray Perform 2DRadon domain transformation and then calculate I using the Grangeat relation xray_Radon , load the trained model and use CT TruePose_Radon ,I xray_Radon Input the trained model, the model predicts I xray The viewpoint pose parameter T coarse , T coarse is the 2D / 3D image registration result based on deep learning prediction, T coarse An example of is as follows.

[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 with a Euclidean distance d from the origin.

[0113] When transforming in the two-dimensional Radon domain, the integral 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 converted into a hyperplane straight line integral,

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

[0121]

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

[0123] Appendix 2: 3DRadon domain transformation formula

[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 with a Euclidean distance d from the origin.

[0127] When transforming in the three-dimensional Radon domain, 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. It represents the angle between the projection of the unit normal vector on the XY plane in the three-dimensional Cartesian coordinate system and the X-axis. (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,

[0134]

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

[0136] Appendix 3: Grangeat relation formula

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

[0138] Assume that the direction of the integral line on the plane detector is the u-axis direction, and the direction perpendicular to it is the v-axis direction.

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

[0140]

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

[0142] All cone beam projection data are multiplied by a cosine function before calculating the data integral. Substituting into the above formula, the data integration on the u-axis is actually:

[0143]

[0144] Substituting variable u for variable θ, we can get formula 4 from formula 3:

[0145]

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

[0147]

[0148] There is also the following relationship:

[0149]

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

[0151] This gives:

[0152]

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

[0154]

[0155] Finally, we can get the relationship between the three-dimensional Radon domain transform and the two-dimensional Radon domain transform of volume data:

[0156]

Claims

1. A 2D / 3D image registration method based on Radon domain transform, characterized in that: The following steps are involved: Step 1, obtaining three-dimensional CT volume data and two-dimensional X-ray images of the target area of ​​the object; Step 2: Preprocess the three-dimensional CT volume data and the two-dimensional X-ray image data respectively; wherein: The preprocessing of the 3D CT volume data is as follows: the image coordinate system of the 3D CT volume data is converted to the 3D coordinate system and the AP view projected from the positive direction of the Y axis to the negative direction of the Y axis is consistent with the direction of the 2D X-ray image, so as to obtain the processed correct CT position image data CT TruePose ; The preprocessing of the two-dimensional X-ray image is as follows: 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 ​​through a normal distribution with a given mean and standard deviation to obtain the preprocessed two-dimensional X-ray image I xray ; Step 3: obtaining a number of existing CT volume data of the same target area of ​​the same type of object as the object in step (1) in the existing CT data set as a training set, and preprocessing the existing CT volume data, including: Step 3.1: Obtain CT volume data - index group CT indexed ; Step 3.2: Sample the true viewpoint pose parameter T from a normal distribution with a given mean and standard deviation gt , through the true viewpoint pose parameter T gt , using the ray casting algorithm to generate DRR image I gt , integrating all DRR images I gt and the corresponding true viewpoint pose parameter T gt Get DRR dataset I drr_dataset ; Step 3.3: DRR dataset I obtained in step 3.2 drr_dataset DRR image I gt Preprocessing: DRR image I gt By normalizing the maximum and minimum pixel values ​​of each image, the resolution of the image is adjusted to the specified size, and the DRR image I gt The pixel values ​​are standardized by a normal distribution with a given mean and standard deviation to obtain the preprocessed DRR dataset I drr_dataset_pre ; Step 4: Perform 3DRadon domain transformation on the existing CT volume data in step 3, and then calculate the 3D training data set by Grangeat relation, i.e., the data set Radon 3D ; Dataset Radon 3D Each data is Radon 3D_input ; For the DRR image I preprocessed in step 3.3 gt Perform 2DRadon domain transformation, and then calculate the 2D training data set by Grangeat relation, that is, the data set Radon 2D ; Dataset Radon 2D Includes Radon 2D_input And its corresponding pose parameter T gt ; Step 5: Build a deep neural network model for 2D / 3D image registration and the loss function used in network training. 2D Randomly select data from Radon 2D_input , and from the Radon dataset 3D Find Radon 2D_input Corresponding Radon 3D_input , Radon 2D_input , Radon 3D_input As model input, the model takes Radon 2D_input , Radon 3D_input Encoded as the corresponding viewpoint pose parameters, the model output is the predicted viewpoint pose parameter T pred ; Step 6: The correct CT position image data CT in step 2 TruePose Perform 3DRadon domain transformation and then calculate CT using the Grangeat relation TruePose_Radon , for the real X-ray image I in step 2 xray Perform 2DRadon domain transformation and then calculate I using the Grangeat relation xray_Radon , load the trained model and use CT TruePose_Radon ,I xray_Radon Input the trained model, the model predicts I xray The viewpoint pose parameter T coarse , T coarse 2D / 3D image registration results based on deep learning prediction.

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, that is, the existing CT volume data is preprocessed and divided into a region with CT value < a HU, a region with a HU < CT value < bHU, and a region with CT value > b HU according to the CT value of each voxel, where HU is the unit of CT value, a and b are both preset values, and an index of the region corresponding to each voxel is made and saved to obtain the preprocessed CT volume data - index group CT indexed .

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 CT generated in step 3.1 is used for each forward projection. indexed The voxel values ​​in the region with CT values ​​> b HU in CT are uniformly multiplied by the bone enhancement coefficient μ, which is randomly sampled 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 physical length D of the diagonal line. CTdiagonal , determine the sampling plane size D CTdiagonal *D CTdiagonal ; According to the physical length of the diagonal D CTdiagonal Determine the discrete sampling spacing of the Radon integration plane, calculate the discrete sampling angle, calculate the sampling coordinates in the physical space of each integration plane, and then convert the sampling coordinates to CT index coordinates, exclude the sampling coordinate points that are out of range, and then accumulate the samples in the CT data to get the integration results, integrate all the sampling points to get the 3DRadon domain transformation results, process them again according to the Grangeat relationship, and standardize them through the normal distribution with a given mean and standard deviation to get the preprocessed 3DRadon data set Radon 3D ; to I drr_dataset_pre Each DRR image I gt The following processing is performed: Calculate the DRR image I gt The real physical size of the image and calculate the DRR image I gt The actual diagonal length is taken as the sampling length D DRRdiagonal ; According to the sampling length D DRRdiagonal Determine the discrete sampling spacing of the Radon integral plane and calculate the discrete sampling angle; calculate the sampling coordinates in the physical space of each integral plane, and then convert the sampling coordinates into DRR image index coordinates, exclude the sampling coordinate points that are out of range, and then accumulate the samples in the DRR image to obtain the integral result, integrate all the sampling points to obtain the 2DRadon domain transformation data, process it again according to the Grangeat relationship, and standardize it through a normal distribution with a given mean and standard deviation to obtain the preprocessed 2DRadon data set Radon 2D .

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 submodules, which are connected in series. The first module consists of a CBR layer and a 3*3*3 maxpool layer. The next four modules have the same minimum unit, which differs only in the number of minimum units and feature extraction depth of each module. The minimum unit consists of two CBR layers plus a Conv layer and a BatchNorm layer and residual connection. The last module is a fully connected layer, which consists of an AvgPool layer and a Linear layer. Its function is to regress the viewpoint pose parameter T from the high-dimensional features extracted from the first five submodules. pted .

6. The 2D / 3D image registration method based on Radon domain transform according to claim 5, characterized in that: There are two loss functions used when building network training. 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 shift MSE loss where dx gt For the current Radon 2D_input The corresponding CT pose is offset from the X axis of the center point, dx pred The X-axis offset predicted by the model; dy gt For the current Radon 2D_input The corresponding CT pose is offset from the Y axis of the center point, dy pred The Y-axis offset predicted by the model; dz gt For the current Radon 2D_input The corresponding CT pose offset relative to the Z axis of the center point, dz pred The Z-axis offset predicted by the model; 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 Angle offset MSE loss 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 yaw angle difference predicted by the model; pitch gt For the current Radon 2D_input The difference between the corresponding CT posture and the initial CT pitch angle, pitch pred is the pitch angle difference predicted by the model; roll gt For the current Radon 2D_input The difference between the corresponding CT posture and the initial CT roll angle, roll pred is the roll angle difference predicted by the model.

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

8. The 2D / 3D image registration method based on Radon domain transform according to claim 5, characterized in that: In the step 5.3, the model is back-propagated and trained by the loss function constructed in step 5.2 to update the model parameters. The 2D / 3D image registration prediction model is trained for a total of num_epoch epochs, each epoch contains num_itr iterations, and each iteration model performs a back-propagation update parameter, and the model file is saved after the training is completed.

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 the 2D / 3D image registration results through a deep neural network model. It is not used as a diagnostic standard and is not for diagnostic or therapeutic purposes.

10. A 2D / 3D image registration system based on Radon domain transform, characterized in that: The following steps are involved: A data acquisition module, used to acquire three-dimensional CT volume data and two-dimensional X-ray images of the target area of ​​the object; The preprocessing module 1 is used to preprocess the three-dimensional CT volume data and the two-dimensional X-ray image data respectively; wherein, The preprocessing of the 3D CT volume data is as follows: the image coordinate system of the 3D CT volume data is converted to the 3D coordinate system and the AP view projected from the positive direction of the Y axis to the negative direction of the Y axis is consistent with the direction of the 2D X-ray image, so as to obtain the processed correct CT position image data CT TruePose ; The preprocessing of the two-dimensional X-ray image is as follows: 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 ​​through a normal distribution with a given mean and standard deviation to obtain the preprocessed two-dimensional X-ray image I xray ; The second preprocessing module is used to obtain a number of existing CT volume data of the same target area of ​​the same type of object as the object in step 1 in the existing CT data set as a training set, and preprocess the existing CT volume data, including: 3.1 Obtaining CT Volume Data - Index Group CT indexed ; 3.2 Sampling the true viewpoint pose parameter T from a normal distribution with a given mean and standard deviation gt , through the true viewpoint pose parameter T gt , using the ray casting algorithm to generate DRR image I gt , integrating all DRR images I gt and the corresponding true viewpoint pose parameter T gt Get DRR dataset I drr_dataset ; 3.3 Acquired DRR Dataset I drr_dataset DRR image I gt Preprocessing: DRR image I gt By normalizing the maximum and minimum pixel values ​​of each image, the resolution of the image is adjusted to the specified size, and the DRR image I gt The pixel values ​​are standardized by a normal distribution with a given mean and standard deviation to obtain the preprocessed DRR dataset I drr_dataset_pre ; The domain transformation processing module is used to perform 3DRadon domain transformation on the existing CT volume data, and then calculate the 3D training data set by Grangeat relationship, that is, the data set Radon 3D ; Dataset Radon 3D Each data is Radon 3D_input ; For the pre-processed DRR image I gt Perform 2DRadon domain transformation, and then calculate the 2D training data set by Grangeat relation, that is, the data set Radon 2D ; Dataset Radon 2D Includes Radon 2D_input And its corresponding pose parameter T gt ; A deep neural network model for 2D / 3D image registration and a loss function module used in network training to extract the Radon dataset 2D Randomly select data from Radon 2D_input , and from the Radon dataset 3D Find Radon 2D_input Corresponding Radon 3D_input , Radon 2D_input , Radon 3D_input As model input, the model takes Radon 2D_input , Radon 3D_input Encoded as the corresponding viewpoint pose parameters, the model output is the predicted viewpoint pose parameter T pred ; Model training and output module, used to convert the correct CT position image data CT TruePose Perform 3DRadon domain transformation and then calculate CT using the Grangeat relation TruePose_Radon , for the real X-ray image I xray Perform 2DRadon domain transformation and then calculate I using the Grangeat relation xray_Radon , load the trained model and use CT TruePose_Radon ,I xray_Radon Input the trained model, the model predicts I xray The viewpoint pose parameter T coarse , T coarse 2D / 3D image registration results based on deep learning prediction.

Images