A 2D / 3D image registration method for robot-assisted surgery system

Abstract

The present invention discloses a 2D / 3D image registration method for a robot-assisted surgery system, which searches for the spatial position relationship between preoperative 3D images (such as CT and MRI images) and intraoperative 2D images (such as X-ray images). The present invention divides 2D / 3D registration into two stages: coarse registration based on deep learning and fine registration based on iteration. The present invention proposes a hyperdimensional rotation SO(4) group isomorphism embedding loss to solve the problems of poor optimization effect of existing deep learning methods on high-dimensional space and low registration accuracy; a new hybrid NCC image similarity metric is proposed, which can increase the capture range of 2D / 3D registration, reduce the probability of falling into local optimality during the registration process, and improve the accuracy of registration based on iterative methods; a sparse rendering method is proposed, which can improve the iteration speed without affecting the registration accuracy. The present invention uses the above technology to improve the accuracy of 2D / 3D registration in a robot-assisted surgery system.

Description

Technical Field

[0001] The present invention belongs to the field of image application, and specifically relates to a multi-stage 2D / 3D registration method for processing two-dimensional X-ray images taken during surgery and three-dimensional CT or MRI volume data taken before surgery in a robot-assisted surgery system. Background Art

[0002] 2D / 3D registration technology is a core component of the robot-assisted surgery system. It is a technology that estimates the relative spatial position relationship between a known 3D structure and its corresponding 2D image. The 3D structure includes the preoperative three-dimensional computed tomography (CT) or magnetic resonance imaging (MRI) image of the patient's skeletal anatomy, and the 2D image includes the intraoperative X-ray fluoroscopic image of the same 3D structure. At present, the existing 2D / 3D registration technology can be divided into iterative optimization-based methods and deep learning-based methods.

[0003] The principle of the 2D / 3D registration method based on iterative optimization is to use the ray casting algorithm to project the 2D digitally reconstructed radiograph (DRR) from the 3D volume data, take the 2D X-ray image as the reference image and the DRR as the floating image, use the search optimization algorithm to iteratively update the DRR from different viewpoints, search for the DRR image with the best similarity measure with the X-ray, and use the viewpoint pose matrix corresponding to the DRR image as the optimization result of 2D / 3D registration. However, the 2D / 3D registration based on iterative optimization is subject to many limitations when applied to clinical practice, including instability in optimization objectives (similarity measure loss function selection), hyperparameter selection and initialization, difficulty in handling inconsistent or multiple objects, and limited single-view X-ray registration performance. Therefore, considering the complexity, variability and high dimensionality of the 2D / 3D registration problem, an important reason for the existence of these problems is that the traditional analytical solution based only on iterative optimization is insufficient in solving such a complex high-dimensional optimization problem.

[0004] At present, the 2D / 3D registration method based on deep learning directly predicts the pose matrix corresponding to the X-ray image through the 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 registration accuracy of the existing deep learning-based method (greater than 10 mm) is much lower than the accuracy of the iterative optimization method (about 2 mm), which is difficult to meet clinical needs. The reason for the low accuracy of the deep learning-based registration method is that the L2 loss is currently mainly used as the distance formula between the true value and the model predicted pose matrix. This loss has a poor optimization effect on the high-dimensional space where the model is located, and it is difficult to better encode the pose of the image through the model. Summary of the invention

[0005] Purpose of the invention: In view of the above shortcomings, the present invention provides a 2D / 3D image registration method for a robot-assisted surgery system, with the aim of providing a 2D / 3D registration algorithm with faster registration speed and higher registration accuracy compared to the prior art.

[0006] Technical solution: To solve the above problems, the present invention discloses a 2D / 3D image registration method for a robot-assisted surgery system, comprising the following steps:

[0007] (1) acquiring three-dimensional CT or MRI volume data and two-dimensional X-ray image data of the target area of ​​the object;

[0008] (2) Preprocessing the three-dimensional volume data and the two-dimensional X-ray image data, dividing the three-dimensional volume data into several regions according to the intensity value interval of each voxel, making and saving the index of the region corresponding to each voxel, and obtaining the preprocessed volume data-index group V indexed ; Sample the two-dimensional X-ray image to a specified size to obtain a preprocessed two-dimensional X-ray image I xray ;

[0009] (3) Volume data-index group V indexed Perform forward projection to generate the DRR dataset I for deep learning training drr_dataset_pre ;

[0010] (4) Construct a deep neural network model for the 2D / 3D coarse registration stage and construct a loss function for network training. drr_dataset_pre Train the deep neural network model; the input of the deep neural network model is the DRR image I drr_input , the deep neural network model transforms the DRR image I drr_input Encoded as the DRR image I drr_input The corresponding viewpoint pose parameters, the output of the deep neural network model is the predicted viewpoint pose parameter T pred ;

[0011] (5) Using 2D X-ray images I xray Input the trained model to predict its position parameters, using the real X-ray image I xray Input the trained deep neural network model, and the deep neural network model predicts I xray The viewpoint pose parameter T coarse As a prediction result;

[0012] (6) The prediction result T coarse Set as the starting position for iterative fine registration to the volume data-index group V indexed Perform forward projection, for the sparse DRR image I drr , after rendering, several DRR local image blocks I are obtained drr|Spars ;

[0013] Construct the iterative optimizer Adam, set the number of iterations to num_itr_refine times, and use the formula in each iteration step:

[0014]

[0015] Calculate the sparse DRR local image block I currently being rendered drr|Spa With the DRR local image block I drr|Spa In the two-dimensional X-ray image I xray The corresponding X-ray local image block I xray|Spars NCC hybrid value; where Sobel(I)=[g x *I,g y *I];

[0016] * represents the convolution operation, g x and g y Represents a 3×3 Sobel operator convolution kernel, which is used to extract the gradient information of the image in the x-axis and y-axis directions. α is a hyperparameter, and i represents the i-th local image block;

[0017] The optimizer Adam updates the viewpoint pose parameters T predicted by the model through gradients pred , projection pose parameter T pred Volume data-index group V indexed Perform forward projection to obtain the sparse DRR image I of the next iteration step drr|Spars When the iteration loop is completed, obtain and save NCC hybrid The projection pose parameter T in the iteration with the largest value pred as the final registration result.

[0018] Furthermore, in step (2), the acquired two-dimensional X-ray image is preprocessed, and the preprocessing is performed in the following order: the X-rays are center-cropped to a uniform size, the X-rays are Gaussian blurred, the image pixel values ​​are logarithmically transformed, the image pixel values ​​are normalized to the range of [0, 1], and the X-ray image is trilinearly interpolated and downsampled to a specified size to obtain a preprocessed two-dimensional X-ray image I xray .

[0019] Furthermore, in step (3), a ray casting algorithm is used to 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 , the DRR image I is generated using the ray casting algorithm through the pose parameters gt ; A total of num_drr_dataset times are sampled to generate num_drr_dataset DRR images, each DRR image I gt With the corresponding true viewpoint pose parameter T gt As labels, a total of num_drr_datasets groups of DRR image-viewpoint pose parameter label pairs are used as training sets for deep learning network models. drr_dataset .

[0020] Furthermore, in each forward projection, the volume data-index group V indexed The voxel values ​​of the bone region in the volume data are uniformly multiplied by the bone enhancement coefficient μ, which is randomly sampled from a uniform distribution within a given range.

[0021] Furthermore, the DRR dataset I obtained in step (3) drr_dataset The preprocessing is performed in the following order: the DRR image is normalized by the maximum and minimum pixel values ​​of each image, the resolution of the image is adjusted to the specified size, and the pixel value of the DRR image is standardized by a normal distribution with a given mean and standard deviation to obtain the preprocessed DRR data set I drr_dataset_pre .

[0022] Furthermore, in step (4), the deep neural network model includes five submodules, and the submodules are connected in a series structure. Each of the first four submodules consists of a downsampling layer and a feature extraction layer. The downsampling layer consists of a strided convolution layer and a Layer-Normalization normalization layer. The feature extraction layer consists of a depthwise separable convolution layer with a convolution kernel size of 7×7 and a reverse bottleneck layer. The reverse bottleneck layer includes two fully connected layers and a GELU activation function. The first four submodules are responsible for extracting high-dimensional features of the DRR images in the training set; the fifth submodule is a fully connected layer, which consists of two fully connected layers, and its function is to regress the viewpoint pose parameter T from the high-dimensional features extracted by the first four submodules. pred .

[0023] Furthermore, in step (4), the deep neural network model is trained by back propagation through constructing a loss function to update the parameters of the deep neural network model. The deep neural network model is trained for num_epoch epochs in total, each epoch contains num_itr iterations, and the deep neural network model performs back propagation to update parameters once in each iteration. After the training is completed, the deep neural network model file is saved.

[0024] Furthermore, in step (4), three loss functions are used to construct network training:

[0025] The first loss function is the input DRR image I drr_input And the viewpoint pose parameter T obtained by forward inference of the model pred By V indexed Projected DRR image I drr_pred The normalized correlation coefficient NCC between them is as follows:

[0026]

[0027] Where Z drr_input and Z drr_pred By I drr_input and I drr_pred The image after normalizing the pixel mean and standard deviation of a single image is as follows:

[0028]

[0029] Wherein, μ(I) represents the average grayscale value of all pixels in image I, σ(I) represents the standard deviation of the grayscale values ​​of all pixels in image I, and Z represents the image after the standardization of image I.

[0030] The second loss function is the input DRR image I drr_input The label truth viewpoint pose parameter T gtAnd the viewpoint pose parameters T predicted by the model pred The dual geodesic distance loss between Where f is the scaling factor constant;

[0031] T gt =[R gt |t gt ]

[0032] T pred =[R pred |t pred ]

[0033]

[0034]

[0035]

[0036] d t (t gt ,t pred )=∥t gt -t pred ∥

[0037] The third loss function is the input DRR image I drr_input The true viewpoint pose parameter T gt And the viewpoint pose parameters T predicted by the model pred The SO(4) group isomorphism between them is embedded in the dual invariant distance loss The formula is as follows:

[0038]

[0039]

[0040]

[0041] Among them, T gt By T gt|SO(4) The formula is mapped to the SO(4) group element T gt|SO(4) , T pred By T pred|SO(4) The formula is mapped to the SO(4) group element T pred|SO(4) , The formula represents the calculation of T gt|SO(4) and T pred|SO(4) Frobenius norm of the difference of matrices.

[0042] Furthermore, in step (6), during the sparse DRR rendering process, for a sparse DRR image I to be rendered with a resolution of 224×224 drr, only 100 DRR local image patches with a resolution of 13×13 are randomly rendered. drr|Spa , a total of 16900 pixels are rendered, and the position index of the center pixel of these 100 local image blocks is retained. center , each local image block is denoted as I patch_drr|Sparse ; For X-ray image I xray , through index center Get I xray The corresponding DRR local image block I drr|Spars 100 X-ray local image patches with a resolution of 13×13 at the position xray|spa , each X-ray local image block is denoted as I patch_xray|sparse .

[0043] Furthermore, in step (2), the three-dimensional volume data is divided into air area, soft tissue area and bone area according to the intensity value of each voxel; the intensity value of the air area is <a; the intensity value of the soft tissue area is: a<intensity value<b; the intensity value of the bone area is >b; where a and b are preset thresholds, and an index of the area corresponding to each voxel is created and saved.

[0044] Beneficial effect: Compared with the prior art, the registration method provided by the present invention has the following advantages: it combines the iterative and deep learning methods into a multi-stage registration method, combining the advantages of high precision of the iterative method and fast running speed of the deep learning method.

[0045] Furthermore, in the registration method, the present invention proposes a new type of loss, maps the SE(3) group pose elements to the SO(4) group, and calculates the distance between the two poses in the SO(4) group, thereby improving the registration accuracy; and, the present invention proposes a new type of hybrid loss function, which mixes and averages the two loss functions of the global NCC with excellent convexity and the patch-GradNCC with significant peaks, and can well represent the cross-modal structural features between the X-ray image and the DRR image while satisfying the excellent convexity and significant optimal peaks, thereby improving the registration accuracy, increasing the 2D / 3D registration capture range, and reducing the probability of the registration falling into the local optimum; and, the present invention uses a sparse rendering method for iteration, rendering only 100 local image blocks of 13×13 pixels each time, and only needs to render 100×13×=16900 pixels at most (assuming that all image blocks do not overlap). Compared with 50176 pixels of the entire image, the rendering amount is only 33.68% of the original image, which greatly improves the iteration speed. Through these technical means, the technical effect of improving the registration speed and improving the registration accuracy can be achieved. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0047] Figure 2 (a) is 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 It is a schematic diagram of the 2D / 3D coarse registration deep neural network model proposed in the technical solution of the present invention.

[0052] Figure 4 (a) is the preprocessed DRR image.

[0053] Figure 4 (b) Yes Figure 4 (a) The corresponding sparsely rendered DRR image.

[0054] Figure 4 (c) is the corresponding sparsely rendered X-ray image at the same iteration. DETAILED DESCRIPTION

[0055] The technical solution of the present invention is further described below in conjunction with the accompanying drawings.

[0056] In order to make the purpose, technical solution and advantages of the present invention more clear, the present invention is further described in detail below in conjunction with the human pelvic structure object CT 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.

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

[0058] Step 1: Obtain 3D CT volume data and 2D X-ray image data of a pelvic object. An example of 3D CT volume data is as follows: Figure 2 (a) shows an example of a two-dimensional X-ray image. Figure 2 (b)

[0059] Step 2, preprocessing the CT data and X-ray data of the pelvic object respectively;

[0060] Step 2.1: Preprocess the 3D CT volume data. Segment the CT volume data into three regions: air (CT value < -800HU), soft tissue (-800HU < CT value < 350HU), and bone (CT value > 350HU) according to the CT value of each voxel (unit: HU). Create and save the index of the region corresponding to each voxel to obtain the preprocessed volume data - index group V. indexed .

[0061] Step 2.2, preprocess the acquired two-dimensional X-ray image. The preprocessing is performed in the following order: crop the center of the X-ray image to a uniform size (resolution 1436×1436 pixels), perform Gaussian blur processing on the X-ray (the size of the Gaussian kernel used is 5×5, and the standard deviation is 1.0), perform logarithmic transformation on the image pixel values, normalize the image pixel values ​​to the range of [0,1] through the maximum grayscale value and the minimum grayscale value of each image, and perform trilinear interpolation downsampling on the X-ray image to the specified size (resolution 224×224 pixels). Finally, the preprocessed two-dimensional X-ray image I is obtained. xray , I xray An image effect like Figure 2 (c) as shown.

[0062] Step 3: Use Siddon's Ray-tracing Algorithm to traverse V indexed Perform forward projection to generate a DRR dataset for deep learning training and preprocess the DRR dataset;

[0063] Step 3.1: Use the ray casting algorithm to indexed Perform forward projection to generate the DRR image. First, sample the true viewpoint pose parameter Y 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(a,b) represents a normal distribution with a mean of a and a standard deviation of b:

[0064] r1=torch.distributions.Normal(0,0.2).sample(), r2=torch.distributions.Normal(0,0.1).sample(), r3=torch.distributions.Normal(0,0.25).sample() , t1=torch.distributions.Normal(10,70).sample(), t2=torch.distributions.Normal(250,90).sample(), t3=torch.distributions.Normal(5,50).sample().

[0065] Through the pose parameter T gt Generating DRR images using ray casting algorithm I gt . A total of 200,000 samples were collected to generate 200,000 DRR images. gt With the corresponding viewpoint pose parameter T gt As labels, a total of 200,000 sets of DRR image-viewpoint pose parameter label pairs are used as the training set I of the deep learning network model drr_dataset .

[0066] Furthermore, in the step 3.1, the volume data-index group V generated in step 2.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.

[0067] Step 3.2: DRR dataset I obtained in step 3.1 drr_dataset Preprocessing is performed in the following order: the DRR images are normalized by the maximum and minimum pixel values ​​of each image, the resolution of the image is adjusted to the specified size (resolution is 224×224 pixels), and 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.

[0068] Step 4: Build a deep neural network model for the 2D / 3D coarse registration stage, build a loss function for network training, and use the DRR dataset I drr_dataset_pre Train the network;

[0069] Step 4.1: Construct a deep neural network model 2D3D-ConvNeXt for the 2D / 3D rough registration stage. The model structure and loss function used are as follows: Figure 3 As shown. Network acquisition data set I drr_dataset_pre DRR image I drr_input As model input. The model consists of five submodules, and the submodules are in series. Each of the first four submodules consists of a downsampling layer and a feature extraction layer. The downsampling layer consists of a strided convolution layer and a Layer-Normalization normalization layer. The feature extraction layer consists of a depth-separable convolution layer with a convolution kernel size of 7*7 and a reverse bottleneck layer. The reverse bottleneck layer includes two fully connected layers and a GELU activation function. The first four submodules are responsible for extracting high-dimensional features of the DRR images in the training set. The fifth submodule is a fully connected layer, which consists of two fully connected layers. Its function is to regress the viewpoint pose parameters T predicted by the network from the high-dimensional features extracted by the first four submodules. pred .

[0070] Step 4.2: Construct three loss functions used for network training.

[0071] The first loss function is the input DRR image I drr_input And the viewpoint pose parameter T obtained by forward inference of the model pred By CT indexed Projected DRR image I drr_pred The normalized correlation coefficient (NCC) between them is as follows:

[0072]

[0073] Where Z drr_input and Z drr_pred By I drr_input and I drr_pred The image after the pixel mean and standard deviation of a single image are standardized is as follows:

[0074]

[0075] Wherein, μ(I) represents the average grayscale value of all pixels in image I, σ(I) represents the standard deviation of the grayscale value of all pixels in image I, and Z represents the image after the image I is standardized.

[0076] The second loss function is the input DRR image I drr_input The label truth viewpoint pose parameter T gt And the viewpoint pose parameters T predicted by the model pred The dual geodesic distance loss between Here, f is the scaling factor constant, which is set to 510.

[0077] T gt =[R gt |t gt ]

[0078] T pred =[R pred |t pred ]

[0079]

[0080]

[0081]

[0082] d t (t gt ,t pred )=∥t gt -t pred ∥

[0083] The third loss function is the input DRR image I drr_input The true viewpoint pose parameter T gt And the viewpoint pose parameters T predicted by the model pred The SO(4) group isomorphism between them is embedded in the dual invariant distance loss The formula is as follows:

[0084]

[0085]

[0086]

[0087] Among them, T gt By T gt|SO(4) The formula is mapped to the SO(4) group element T gt|SO(4) , T pred By T pred|SO(4) The formula is mapped to the SO(4) group element T pred|SO(4) , The formula represents the calculation of T gt|SO(4) and T pred|SO(4) Frobenius norm of the difference of matrices.

[0088] Step 4.3: Use the prepared DRR dataset I drr_dataset_pre Train the network. The input of the model is the DRR image I drr_input , the model will I drr_input Encoded as its corresponding viewpoint pose parameter, the model output is the predicted viewpoint pose parameter T pred, where T pred It is a rigid pose transformation matrix of size 4×4.

[0089] Furthermore, in step 4.3, Figure 3 The size of the output features of each submodule is shown in , where B represents the batch size used for network training. B is set to 2 during model training. The size of the DRR image input to the model is (2,1,224,224), which means that the model uses two DRR images as input in one iteration. The pose parameter matrix label examples corresponding to the two DRR images are as follows:

[0090]

[0091]

[0092] The feature size after the first submodule is (2,96,56,56), the feature size after the second submodule is (2,192,28,28), the feature size after the third submodule is (2,384,14,14), the feature size after the fourth submodule is (2,768,7,7), and after the fifth submodule, two features of size (2,3) are formed. These two features are concatenated into a feature of (2,6), and then exponentially mapped to a predicted pose parameter matrix T of size (2,4,4). pred , represents two rigid transformation pose parameter matrices of size 4×4, T pred An example is as follows:

[0093]

[0094]

[0095] Furthermore, in the step 4.3, the model is back-propagated and trained using the loss function constructed in step 4.2 to update the model parameters. A three-dimensional CT data prediction model is trained for a total of 1000 epochs, each epoch contains 100 iterations, and the Adam optimizer is used for optimization. The model performs back-propagation and updates the parameters each time it is iterated. After the training is completed, the model ".pth" file is saved.

[0096] Step 5: Using real X-ray images I xray Input the trained model to predict its position parameters, using the real X-ray image I xray Input the trained model, the model predicts I xray The viewpoint pose parameter T coarse , T coarse is the 2D / 3D registration result predicted based on deep learning, T coarse An example is as follows:

[0097]

[0098] Step 6: In the fine registration stage, a gradient-based iterative method is used to perform 2D / 3D registration, and a similarity measurement loss function for registration is constructed. To speed up the registration, a DRR image is generated by a sparse rendering method in each iteration, and iterative optimization is performed through the built-in optimizer Adam in the PyTorch framework to find the optimal registration pose parameters.

[0099] Step 6.1: The deep learning model prediction result T obtained in step 4.2 coarse Set as the starting position of iterative fine registration, through T coarse CT indexed Perform forward projection and render sparse DRR.

[0100] Furthermore, in the sparse rendering process, for a DRR projection image I with a resolution of 224×224 drr , only 100 local image patches with a resolution of 13×13 are randomly rendered I drr|Sparse , a total of 16900 pixels are rendered, and the position index of the center pixel of these 100 local image blocks is retained. center , each local image block is denoted as I patch_drr|Sparse For X-ray image I xray , through index center Get I xray 100 local image blocks with a resolution of 13×13 at the corresponding positions in xray|Sparse , each local image block is denoted as I patch_xray|Sparse A corresponding DRR original image and its sparse rendering image are respectively as follows: Figure 4 (a) Figure 4 (b) as shown.

[0101] Step 6.2: Construct a new hybrid NCC loss function for iterative registration: Hybrid Normalized Cross Correlation (Hybrid-NCC), denoted as NCC hybrid Among them, * represents the convolution operation, g x and g y It represents a 3*3 Sobel operator convolution kernel, which is used to extract the gradient information of the image in the x-axis and y-axis directions. α is a hyperparameter and is set to 0.5.

[0102]

[0103] Sobel(I)=[g x *I,g y*I]

[0104]

[0105] Step 6.3: Build the iterative optimizer Adam, set the number of iterations to 200, and use the above NCC in each iteration step hybrid (I drr|Sparse ,I xray|Sparse ) formula to calculate the currently rendered sparse DRR image I drr|Spars With sparse X-ray image I xray|Sparse NCC hybrid The optimizer Adam updates the projection pose parameter T through the gradient pred , T pred By CT indexed Project the sparse DRR image I of the next iteration step drr|Spa , according to the current I drr|Spa Update the current sparse X-ray image I by the position of the local image block in xray|Spars , the sparsely rendered X-ray image in one iteration step is as follows Figure 4 (c) is shown. The goal of the optimizer is to make the current I drr|Sparse with I xray|Sparse NCC hybrid When the iteration loop is completed, obtain and save NCC hybrid T in the iteration with the largest value pred As the final registration result of the multi-stage medical image 2D / 3D registration method for a robot-assisted surgery system proposed by the present invention. pred An example is as follows:

[0106]

[0107] The above registration algorithm was used for verification on a public pelvic dataset (6 CT data, a total of 366 X-ray data) and a private pelvic dataset (1 CT data, a total of 20 X-ray data). The verification results can achieve an average mTRE (mm) of 0.71±2.87 mm (mTRE refers to the mean target registration error, which calculates the spatial distance between the markers in the dataset and the registered markers).

[0108] After initialization using the deep learning registration method of the present invention, an iterative method is used for fine registration, and finally a 93.45% successful registration rate within 1 mm (generally referred to as a sub-millimeter successful registration rate) can be achieved. The successful registration rate here means that, for example, there are 100 X-ray images to be registered, and after registration, the registration error of 90 images is within 1 mm, then the sub-millimeter successful registration rate is 90%.

[0109] At the same time, the registration method of the present invention can achieve an intraoperative registration time of less than 10 seconds.

[0110] In summary, the present invention aims to address the problem that the optimization time of the existing 2D / 3D registration algorithm based on iterative optimization is too long, while the accuracy of the method based on deep learning cannot meet the requirements. A multi-stage registration strategy is proposed, which divides the 2D / 3D registration into two stages: coarse registration based on deep learning and fine registration based on iterative method. The advantages of the two methods are combined to obtain a 2D / 3D registration algorithm with fast speed and high accuracy.

[0111] In order to solve the problem that the existing image similarity measurement loss function has poor convexity, which leads to the problem that the registration is easy to fall into the local optimum and the registration capture range is small, a new hybrid loss function is proposed. This loss function can well represent the cross-modal structural features between X-ray images and DRR images while meeting the requirements of excellent convexity and significant optimal peak, thereby improving the registration accuracy, increasing the 2D / 3D registration capture range, and reducing the probability of the registration falling into the local optimum.

[0112] In order to solve the problem of poor accuracy of existing deep learning 2D / 3D registration methods, the present invention proposes a novel hyper-dimensional rotation SO(4) group isomorphic embedding loss function to train the network, thereby improving the accuracy of the deep learning-based 2D / 3D registration method.

[0113] In order to solve the problem of long optimization time of current iteration-based 2D / 3D registration methods, the present invention proposes a sparse rendering method. When generating DRR in each iteration, only a part of the pixels of the entire DRR image (no more than 34% of the total number of pixels of a DRR image) is used, which greatly improves the iteration speed while maintaining the accuracy of the 2D / 3D registration results.

Claims

1. A 2D / 3D image registration method for a robot-assisted surgery system, characterized in that: Including the following steps: (1) Obtain three-dimensional CT or MRI volume data and two-dimensional X-ray image data of the object target area; (2) Preprocessing the three-dimensional volume data and the two-dimensional X-ray image data, dividing the three-dimensional volume data into several regions according to the intensity value interval of each voxel, making and saving the index of the region corresponding to each voxel, and obtaining the preprocessed volume data-index group V indexed ; Sample the two-dimensional X-ray image to a specified size to obtain a preprocessed two-dimensional X-ray image I xray ; (3) Volume data-index group V indexed Perform forward projection to generate the DRR dataset I for deep learning training drr_dataset_pre ; In step (3), use the ray casting algorithm to 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 , the DRR image I is generated using the ray casting algorithm through the pose parameters gt ; A total of num_drr_dataset times are sampled to generate num_drr_dataset DRR images, each DRR image I gt With the corresponding true viewpoint pose parameter T gt As labels, a total of num_drr_datasets groups of DRR image-viewpoint pose parameter label pairs are used as training sets for deep learning network models. drr_dataset ; At each forward projection, pass the volume data-index group V indexed The voxel values ​​of the bone region in the volume data are uniformly multiplied by the bone enhancement coefficient β, and μ is randomly sampled from a uniform distribution in a given range; for the acquired DRR dataset I drr_dataset The preprocessing is performed in the following order: the DRR image is normalized by the maximum and minimum pixel values ​​of each image, the resolution of the image is adjusted to the specified size, and the pixel value of the DRR image is standardized by a normal distribution with a given mean and standard deviation to obtain the preprocessed DRR data set I drr_dataset_pre ; (4) Construct a deep neural network model for the 2D / 3D image coarse registration stage and construct a loss function used in network training. drr_dataset_pre Train the deep neural network model; the input of the deep neural network model is the DRR image I drr_input , the deep neural network model transforms the DRR image I drr_input Encoded as the DRR image I drr_input The corresponding viewpoint pose parameters, the output of the deep neural network model is the predicted viewpoint pose parameter T pred ; In step (4), the deep neural network model includes five submodules, and the submodules are connected in a series structure. Each of the first four submodules consists of a downsampling layer and a feature extraction layer. The downsampling layer consists of a strided convolution layer and a Layer-Normalization normalization layer. The feature extraction layer consists of a depthwise separable convolution layer with a convolution kernel size of 7×7 and a reverse bottleneck layer. The reverse bottleneck layer includes two fully connected layers and a GELU activation function. The first four submodules are responsible for extracting high-dimensional features of the DRR images in the training set; The fifth submodule is a fully connected layer, which consists of two fully connected layers. Its function is to regress the viewpoint pose parameters T from the high-dimensional features extracted from the first four submodules. pred ; The deep neural network model is trained by back propagation through constructing loss functions, and the parameters of the deep neural network model are updated. The deep neural network model is trained for num_epoch epochs in total, and each epoch contains num_itr iterations. The deep neural network model performs back propagation to update parameters once in each iteration, and the deep neural network model file is saved after training is completed; Three loss functions are used to construct network training: The first loss function is the input DRR image I drr_input And the viewpoint pose parameter T obtained by forward inference of the model pred By V indexed Projected DRR image I drr_pred The normalized correlation coefficient NCC between them is as follows: Where Z drr_input and Z drr_pred By I drr_input and I drr_pred The image after normalizing the pixel mean and standard deviation of a single image is as follows: Among them, μ(I) represents the average value of the gray levels of all pixels in the calculated image, σ(I) represents the standard deviation of the gray levels of all pixels in the calculated image I, and Z represents the image after normalization; The second loss function is the input DRR image I drr_input The label truth viewpoint pose parameter T gt And the viewpoint pose parameters T predicted by the model pred The dual geodesic distance loss between Where f is the scaling factor constant; T gt =[R gt |t gt ] T pred =[R pre d|t pred ] d t (t gt ,t pred )=||t gt -t pred || The third loss function is the input DRR image I drr_input The true viewpoint pose parameter T gt And the viewpoint pose parameters T predicted by the model pred The SO(4) group isomorphism between them is embedded in the dual invariant distance loss The formula is as follows: Among them, T gt By T gt|SO(4) The formula is mapped to the SO(4) group element T gt|SO(4) , T pred By T pred|SO(4) The formula is mapped to the SO(4) group element T pred|SO(4) , The formula represents the calculation of T gt|SO(4) and T pred|SO(4) Frobenius norm of the difference of matrices; (5) Using 2D X-ray images I xray Input the trained model to predict its position parameters, using the real X-ray image I xray Input the trained deep neural network model, and the deep neural network model predicts I xray The viewpoint pose parameter T coarse As a prediction result; (6) The prediction result T coarse Set as the starting position for iterative fine registration to the volume data-index group V indexed Perform forward projection, for the sparse DRR image I drr , after rendering, several DRR local image blocks I are obtained drr|Sparse ; Construct an iterative optimizer Adam, set the number of iterations to num_itr_refine times, and in each iteration step, use the formula: Calculate the sparse DRR local image block I currently being rendered drr|Sparse With the DRR local image block I drr|Spars In the two-dimensional X-ray image I xray The corresponding X-ray local image block I xray|Spars NCC hybrid value; where Sobel(I)=[g x *I, g y *I]; * represents the convolution operation, g x and g y Represents a 3×3 Sobel operator convolution kernel, which is used to extract the gradient information of the image in the x-axis and y-axis directions. α is a hyperparameter, and i represents the i-th local image block; The optimizer Adam updates the viewpoint pose parameters T predicted by the model through gradients pred , projection pose parameter T pred Volume data-index group V indexed Perform forward projection to obtain the sparse DRR image I of the next iteration step drr|Spars When the iteration loop is completed, obtain and save NCC hybrid The projection pose parameter T in the iteration with the largest value pred as the final registration result.

2. The 2D / 3D image registration method for a robot-assisted surgery system according to claim 1, characterized in that: In step (2), the acquired two-dimensional X-ray image is preprocessed, and the preprocessing is performed in the following order: the X-rays are center-cropped to a uniform size, the X-rays are Gaussian blurred, the image pixel values ​​are logarithmically transformed, the image pixel values ​​are normalized to the range of [0, 1], and the X-ray image is trilinearly interpolated and downsampled to a specified size to obtain a preprocessed two-dimensional X-ray image I xray .

3. The 2D / 3D image registration method for a robot-assisted surgery system according to claim 1, characterized in that: In step (6), during the sparse DRR rendering process, for a sparse DRR image I with a resolution of 224×224 to be rendered, drr , only 100 DRR local image patches with a resolution of 13×13 are randomly rendered. drr|Spa , a total of 16900 pixels are rendered, and the position index of the center pixel of these 100 local image blocks is retained. center , each local image block is denoted as I patch_drr|Spar ; For X-ray image I xray , through index center Get I xray The corresponding DRR local image block I drr|Sparse 100 X-ray local image patches with a resolution of 13×13 at the position xray|Spa , each X-ray local image block is denoted as I patch_xray|Sparse .

4. The 2D / 3D image registration method for a robot-assisted surgery system according to claim 1, characterized in that: In step (2), the volume data is segmented into an air region, a soft tissue region, and a bone region according to the intensity value of each voxel; the intensity value of the air region < a; the intensity value of the soft tissue region is: a < intensity value < b; the intensity value of the bone region > b; where a and b are preset thresholds, and the index corresponding to each voxel region is made and saved.

Images