DTI registration method based on dual-channel parallel network

Abstract

The invention provides a DTI registration method based on a dual-channel parallel network, and relates to the field of medical image processing. A DTI derivative image registration method is combined with a deep learning network, and a DTI derivative scalar FA image and a main feature vector are input into a PD network for training, so that the registration performance is improved while the operand is not increased; a main feature vector in fiber bundle tracking is introduced into a PD network, and an FA image is assisted to predict a more accurate deformation field; for the complexity of the PD network and the increase of input information, an attention mechanism module is introduced into a jump connection part in the middle of each sub-network, so that the PD network suppresses irrelevant regions in extracted features, highlights specific significant features, and improves the registration precision. According to the method, the registration speed and the registration performance are checked through multiple evaluation methods; the PD network provided by the invention has generalization.

Description

Technical Field

[0001] The present invention relates to the field of medical image processing, and in particular to a DTI registration method based on a dual-channel parallel network. Background Art

[0002] At present, brain lesions are still one of the key directions of medical research. Diffusion Tensor Magnetic Resonance Image (DT-MRI, DTI for short) can provide clear brain anatomical structure and rich white matter fiber bundle information, truly reflect the direction of brain anatomical connectivity, and is of great help to doctors in judging the condition after image registration. At present, traditional DTI registration algorithms can be divided into tensor-based registration and DTI-derived image (scalar image)-based registration. Among them, the Diffeomorphic Deformable Tensor Registration (DTI-TK) algorithm has the best registration performance among the tensor-based registration algorithms, and the Symmetric Image Normalization (SyN) algorithm has the best registration performance among the DTI-derived image-based registration algorithms. However, these two types of registration methods have defects such as high computational cost and insufficient use of tensor information, which make them difficult to apply to actual DTI registration.

[0003] The deep learning-based diffusion tensor magnetic resonance image registration method is an efficient registration method. At present, deep learning technology has been applied to the registration of medical images such as magnetic resonance images (MRI) and computed tomography (CT), but there has been no research results on direct registration of DTI data. The traditional registration method based on DTI-derived images only uses fractional anisotropy (FA) images for training to achieve DTI registration. However, the registration method based on DTI-derived images itself has the problem of information loss: the diffusion tensor of DTI contains information on three eigenvalues ​​and three eigenvectors at the same time, while FA only contains information on DTI eigenvalues. Therefore, the traditional method does not use the tensor direction information, and there is room for further improvement in registration performance. In addition, considering that in the process of DTI fiber bundle tracking, the information of the main eigenvector is generally used, while the complementarity of different information of FA and the main eigenvector in characterizing DTI is ignored, the present invention combines them to improve the performance of registration; on the other hand, although the main eigenvector is not a three-dimensional scalar image, the main eigenvector does not need to satisfy the symmetric positive definiteness, so the present invention introduces the main eigenvector into the registration network, and uses it as the network input together with the FA image to assist the FA image in predicting a more accurate deformation field. In order to realize this network input, two VoxelMorph networks that can only input a pair of 3D medical images are improved and connected in parallel as sub-networks to construct a registration method of a dual-channel parallel (Parallel Dual-channel, i.e. PD) network, a registration method of a dual-channel parallel network (DTI Registration method Based on PD), referred to as DTI RBPD. Summary of the invention

[0004] In order to overcome the above technical problems, the present invention proposes a DTI registration method based on a dual-channel parallel network, which combines the registration method of DTI-derived images with a deep learning network, uses the derived scalar FA map and main eigenvector of DTI to input into the dual-channel parallel (Parallel Dual-channel, PD) network for training, and improves the registration performance without increasing the amount of calculation; introduces the main eigenvector in the fiber bundle tracking into the PD network to assist the FA image in predicting a more accurate deformation field; in view of the increase in the complexity of the PD network and the input information, an attention mechanism module is introduced into the jump connection part in the middle of each sub-network, so that the PD network suppresses irrelevant areas in the extracted features, and at the same time highlights specific significant features to improve the registration accuracy. The present invention verifies the registration speed and registration performance through a variety of evaluation methods; the PD network of the present invention has generalization.

[0005] The present invention proposes a DTI registration method based on a dual-channel parallel network, and the specific steps of the technical solution adopted are as follows:

[0006] Step 1, determine the data set: select the IXI (Information eXtraction from Images) data set, preprocess the IXI data set to obtain the DTI data set, and divide the DTI data set into a training set and a test set;

[0007] Step 2: Build a PD network;

[0008] Step 3, training the PD network;

[0009] Step 4, test the PD network;

[0010] Step 5, evaluate the registered DTI;

[0011] Step 6: Verify the generalization of DTI RBPD on different datasets: Use the dHCP (developing Human Connectome Project) dataset and repeat steps 1 to 5 for verification.

[0012] Further, in step 1, the IXI dataset preprocessing method is as follows:

[0013] First, the diffusion image in the IXI dataset is processed by eddy current correction to obtain Image 1; then the skull in Image 1 is removed to obtain Image 2; finally, a DTI dataset with an individual resolution of 128×128×64 is obtained by tensor calculation; the DTI dataset is divided into a training set and a test set in a ratio of 8:1;

[0014] A DTI with excellent morphology and image quality was manually observed and randomly selected from the DTI dataset as the fixed DTI in the registration process, and the remaining DTIs were the DTIs to be registered. The derived FA images and main eigenvectors of each DTI in the DTI dataset were calculated using FSL software. The FA image included the fixed FA and the FA to be registered. The main eigenvector included the fixed main eigenvector and the main eigenvector to be registered. The FA image pair was the fixed FA and the FA to be registered, and the main eigenvector pair was the fixed main eigenvector and the main eigenvector to be registered.

[0015] Further, in step 2, the PD network construction method is as follows:

[0016] The PD network includes two identical sub-networks and a spatial transformation module; the sub-network is a VoxelMorph network that introduces an attention mechanism; the sub-network for inputting FA image pairs is the upper sub-network, and the sub-network for inputting main feature vector pairs is the lower sub-network;

[0017] S2.1, input the FA image pair into the upper sub-network, and input the main feature vector pair of the DTI to be registered and the fixed DTI into the lower sub-network; the FA image and DTI are output as two deformation fields through the two sub-networks respectively;

[0018] S2.2, the two deformation fields pass through a convolution layer and an activation function at the same time, and the output is a fused deformation field, and the number of output channels of the convolution layer is 3;

[0019] S2.3, fusion of the deformation field and the DTI to be registered into the spatial transformation module to obtain the deformed DTI;

[0020] S2.4, the deformed DTI is output as the registered DTI after tensor redirection;

[0021] S2.5, calculate the loss function L of the PD network.

[0022] Further, the sub-network includes an encoder, a tandem skip connection and a decoder;

[0023] The encoder uses a three-dimensional convolutional layer with a convolution kernel size of 3×3×3 and a step size of 2, and the convolutional layer is followed by a LeakyReLU activation function; the encoder reduces the spatial dimension of the feature map through convolution to obtain DTI shallow features;

[0024] The decoder uses a three-dimensional convolutional layer with a convolution kernel size of 3×3×3 and a step size of 2. The convolutional layer is followed by a LeakyReLU activation function. The decoder enlarges the spatial dimension of the feature map through upsampling and convolutional layers to obtain DTI deep features. The spatial dimension of the decoder changes in the opposite direction to that of the encoder. The number of output channels of the decoder is 3, corresponding to the number of channels of the predicted deformation field.

[0025] The serial jump connection is a convolutional block attention module (CBAM) in which channel attention and spatial attention are connected in series. CBAM assigns weights to deep features and shallow features, thereby suppressing irrelevant areas in deep features and shallow features, while highlighting specific significant features, so that deep features and shallow features are fused to achieve the purpose of feature fusion.

[0026] Furthermore, the loss function L is calculated as follows:

[0027] The loss function of the PD network includes the similarity loss L sim and smoothing loss L smooth Two parts; L sim The similarity loss is calculated by fixing the FA image and the deformed FA image, L smoothThe smoothing loss is calculated by the deformation field; the calculation formula of the PD network loss function L is as follows:

[0028]

[0029] Among them, f and m are the three-dimensional space domain Ω∈R n In the image, f represents the fixed DTI, m represents the DTI to be registered, represents the deformation field, represents the deformed image output by the spatial transformation module, and ξ is the regularization parameter.

[0030] Since FA images clearly retain the tissue structure and position information of DTI, traditional registration methods based on DTI-derived images all use FA images to iteratively calculate the deformation field. Therefore, the principal eigenvector is mainly used as complementary information to cooperate with FA images to obtain a more accurate deformation field, and does not participate in the calculation of the loss function.

[0031] Further, in step 3: the training method of the PD network is as follows:

[0032] The FA image of the DTI to be registered and the FA image of the fixed DTI in the training set are stacked as the input of the upper sub-network; the main eigenvector of the DTI to be registered and the main eigenvector of the fixed DTI in the training set are stacked as the input of the lower sub-network; the two sub-networks output the FA predicted deformation field and the main eigenvector predicted deformation field respectively, and the FA predicted deformation field and the main eigenvector predicted deformation field are input into the convolution layer together, and after the LeakyReLU activation function, a fused deformation field with 3 channels is obtained; the fused deformation field and the FA to be registered are output as the deformed FA after passing through the spatial transformation module, which is used to calculate the loss function and then train the network.

[0033] Further, in step 4: the testing method of the PD network is as follows:

[0034] The FA images of the DTI to be registered and the fixed DTI in the test set are input into the upper sub-network of the trained PD network, and the main feature vectors of the DTI to be registered and the fixed DTI in the training set are input into the lower sub-network of the trained PD network; two deformation fields are obtained through the two sub-networks of the PD network, and the two deformation fields pass through a convolution layer and an activation function to output a fused deformation field; the fused deformation field and the DTI to be registered are input into the spatial transformation module to obtain the deformed DTI, and the deformed DTI uses the ANTs toolkit for tensor redirection, and finally the registered DTI is obtained.

[0035] Furthermore, in step 5, the evaluation indexes of the registered DTI include subjective visual evaluation, registration time comparison, region-based evaluation and fiber bundle-based evaluation.

[0036] Further, in step 5, the evaluation method of the evaluation index is as follows:

[0037] The DTI output from step 4 and the fixed DTI selected from the IXI dataset were subjected to subjective visual evaluation, registration time comparison, region-based evaluation, and fiber bundle-based evaluation.

[0038] The region-based evaluation method includes two parts: tensor-based evaluation measure and scalar-based evaluation measure. The tensor-based evaluation measure and the scalar-based evaluation measure are compared through the data average value.

[0039] The tensor-based evaluation measures include tensor Euclidean distance and tensor deviation Euclidean distance; the scalar-based evaluation measures include bidirectional consistency, eigenvalue-eigenvector pair overlap rate, cross-correlation, and the angle between the main eigenvector pairs;

[0040] (I) Tensor Euclidean Distances (ED)

[0041] The smaller the ED value, the better the registration performance and the more aligned the biological tissue structure in DTI. The ED calculation formula is as follows:

[0042]

[0043] Where D1 and D2 represent two diffusion tensors respectively. The diffusion tensor is a 3×3 symmetric positive definite matrix. Tr() represents the trace of the tensor. The calculation formula of Tr() is as shown in formula (3):

[0044] Tr(D)=D xx +D yy +D zz (3)

[0045] Among them, D xx , D yy and D zz They are the three elements on the diagonal of the diffusion tensor D;

[0046] (ii) Euclidean distance of the tensor deviation (DD)

[0047] The smaller the DD value, the better the registration performance and the more aligned the biological tissue structure in DTI. The DD calculation formula is as follows:

[0048]

[0049] Among them, ||·|| c Indicates calculation of Euclidean distance;

[0050] (III) Dyadic Coherence (κ)

[0051] The registration results are evaluated using the bidirectional consistency of the tensor. The changes in the biological tissue structure are evaluated by measuring the changes in the direction of the DTI main eigenvector. The κ calculation formula is as follows:

[0052]

[0053] Among them, β i It is called a dyadic tensor; e1, e2, and e3 are the three eigenvectors of the diffusion tensor D; e1 is the main eigenvector; e 1x ,e 1y ,e 1z are the three components of the principal eigenvector respectively; is the principal eigenvector of the i-th subject; the larger the κ value, the better the performance of the registration algorithm;

[0054] (IV) Overlap of Eigenvalue-Eigenvector Pairs (OVL)

[0055] The larger the OVL value, the better the performance of the registration algorithm. The calculation formula of OVL is as follows:

[0056]

[0057] in, and are the jth eigenvalue and eigenvector of the i-th voxel of the registered DTI, respectively; and are the j-th eigenvalue and eigenvector of the i-th voxel of the fixed DTI, respectively;

[0058] 5. Cross Correlations (CC)

[0059] The cross-correlation between voxels is calculated using the scalar derivative measure of DTI. The larger the cross-correlation value, the higher the similarity of the images. The calculation formula of the cross-correlation CC is as follows:

[0060]

[0061] Where v represents the voxel of all images; X1(v) and X2(v) are scalar images; the value range of CC is [0,1];

[0062] The correlation between FA and Tr of the registered DTI and the fixed DTI was calculated respectively; CC Tr Calculation method and CCFA Same; FA's correlation CC FA The calculation formula is as follows (10):

[0063]

[0064] (VI) Angle between the principal eigenvector pairs

[0065] The angle (α) of the main eigenvector is:

[0066]

[0067] Among them, e′1 represents the main eigenvector of the fixed DTI, e1 represents the main eigenvector of the registered DTI, and the smaller the α value, the better the performance of the registration algorithm.

[0068] (VII) Evaluation based on fiber bundles: The fiber bundles of the registered DTI image are obtained by mapping the fiber bundles of the fixed image, and then the characteristic average is calculated and the average characteristic curve along the fiber bundle is plotted; the larger the normative correlation of the FA curves before and after registration, the better; the more similar the means are, the better; and the smaller the variance is, the better; the FA characteristic curves of the selected fiber bundles include the genu of the corpus callosum (Genu), the splenium of the corpus callosum (Splenium), the left-superior thalamic radiation (L-str), the right-superior thalamic radiation (R-str), the left corticospinal / corticobulbar tracts (L-cst), and the right corticospinal / corticobulbar tracts (R-cst).

[0069] Technical effects of the present invention:

[0070] (1) The present invention combines the registration method of DTI-derived images with a deep learning network, and uses the DTI-derived scalar FA map and the main eigenvector to input into the PD network for training, thereby improving the registration performance without increasing the amount of computation;

[0071] (2) The present invention introduces the main eigenvector in fiber bundle tracking into the PD network to assist FA images in predicting more accurate deformation fields;

[0072] (3) In view of the increase in PD network complexity and input information, the present invention introduces an attention mechanism module in the jump connection part in the middle of each sub-network, so that the PD network suppresses irrelevant areas in the extracted features and highlights specific significant features, thereby improving the registration accuracy;

[0073] (4) The present invention tests the registration speed and registration performance through a variety of evaluation methods, and the PD network of the present invention has generalization. BRIEF DESCRIPTION OF THE DRAWINGS

[0074] Figure 1 It is the overall flow chart of the present invention;

[0075] Figure 2 It is an axial section diagram of the data preprocessing result of the present invention;

[0076] Figure 2 In the figure, (a) is the diffusion image; (b) is the image 1 after eddy current correction; (c) is the image 2 after skull removal; (d) is the DTI image;

[0077] Figure 3 This is a dual-channel parallel network structure diagram of the present invention;

[0078] Figure 4 It is a subnetwork structure diagram of the dual-channel parallel network of the present invention;

[0079] Figure 5 This is the structural diagram of the convolutional attention module of the present invention;

[0080] Figure 6 This is a visual comparison diagram of the axial position of the registration result of the present invention;

[0081] Figure 6 In the figure, (a) is the fixed DTI under axial vision; (b) is the DTI to be registered under axial vision; (c) is the registration result of the SyN method under axial vision; (d) is the registration result of the DTI-TK method under axial vision; (e) is the registration result of the DTI RBV method under axial vision; (f) is the registration result of the DTI RBPD method under axial vision;

[0082] Figure 7 This is a coronal visual comparison diagram of the registration result of the present invention;

[0083] Figure 7 In the figure, (a) is the fixed DTI under coronal vision; (b) is the DTI to be registered under coronal vision; (c) is the SyN registration result under coronal vision; (d) is the DTI-TK registration result under coronal vision; (e) is the DTI RBV registration result under coronal vision; (f) is the DTI RBPD registration result under coronal vision;

[0084] Figure 8 This is a sagittal visual comparison diagram of the registration result of the present invention;

[0085] Figure 8 In the figure, (a) is the fixed DTI under sagittal vision; (b) is the DTI to be registered under sagittal vision; (c) is the SyN registration result under sagittal vision; (d) is the DTI-TK registration result under sagittal vision; (e) is the DTI RBV registration result under sagittal vision; (f) is the DTI RBPD registration result under sagittal vision; DETAILED DESCRIPTION

[0086] The present invention proposes a DTI registration method based on a dual-channel parallel network, and the specific steps of the technical solution adopted are as follows:

[0087] Step 1, determine the data set: select the IXI (Information eXtraction from Images) data set, preprocess the IXI data set to obtain the DTI data set, and divide the DTI data set into a training set and a test set;

[0088] Step 2: Build a PD network;

[0089] Step 3, training the PD network;

[0090] Step 4, test the PD network;

[0091] Step 5, evaluate the registered DTI;

[0092] Step 6: Verify the generalizability of DTI RBPD on different datasets: Use the dHCP dataset and repeat steps 1 to 5 for verification. Figure 1 shown.

[0093] Further, in step 1, the IXI dataset preprocessing method is as follows:

[0094] First, the diffusion image in the IXI dataset is processed by eddy current correction to obtain image 1; then the skull in image 1 is removed to obtain image 2; finally, a DTI dataset with an individual resolution of 128×128×64 is obtained by tensor calculation; Figure 2 As shown in Figure 2, the DTI dataset is divided into a training set and a test set in a ratio of 8:1.

[0095] A DTI with excellent morphology and image quality was manually observed and randomly selected from the DTI dataset as the fixed DTI in the registration process, and the remaining DTIs were the DTIs to be registered. The derived FA images and main eigenvectors of each DTI in the DTI dataset were calculated using FSL software. The FA image included the fixed FA and the FA to be registered. The main eigenvector included the fixed main eigenvector and the main eigenvector to be registered. The FA image pair was the fixed FA and the FA to be registered, and the main eigenvector pair was the fixed main eigenvector and the main eigenvector to be registered.

[0096] Further, in step 2, the PD network construction method is as follows:

[0097] The PD network includes two identical sub-networks and a spatial transformation module; the sub-network is a VoxelMorph network that introduces an attention mechanism; the sub-network for inputting FA image pairs is the upper sub-network, and the sub-network for inputting main feature vector pairs is the lower sub-network;

[0098] S2.1, input the FA image pair into the upper sub-network, and input the main feature vector pair of the DTI to be registered and the fixed DTI into the lower sub-network; the DTI to be registered and the fixed DTI respectively output two deformation fields through the two sub-networks;

[0099] S2.2, the two deformation fields pass through a convolution layer and an activation function, and the output is a fused deformation field, and the number of output channels of the convolution layer is 3;

[0100] S2.3, fusion of the deformation field and the DTI to be registered into the spatial transformation module to obtain the deformed DTI;

[0101] S2.4, the deformed DTI is output as the registered DTI after tensor redirection;

[0102] S2.5, calculate the loss function L of the PD network. Figure 3 shown.

[0103] Further, the sub-network includes an encoder, a tandem skip connection and a decoder;

[0104] The encoder uses a three-dimensional convolutional layer with a convolution kernel size of 3×3×3 and a step size of 2, and the convolutional layer is followed by a LeakyReLU activation function; the encoder reduces the spatial dimension of the feature map through convolution to obtain DTI shallow features;

[0105] The decoder uses a three-dimensional convolutional layer with a convolution kernel size of 3×3×3 and a step size of 2. The convolutional layer is followed by a LeakyReLU activation function. The decoder enlarges the spatial dimension of the feature map through upsampling and convolutional layers to obtain DTI deep features. The spatial dimension of the decoder changes in the opposite direction to that of the encoder. The number of output channels of the decoder is 3, corresponding to the number of channels of the predicted deformation field.

[0106] The series jump connection is a convolutional block attention module (CBAM) that connects channel attention and spatial attention in series. CBAM assigns weights to deep features and shallow features, thereby suppressing irrelevant areas in deep features and shallow features, while highlighting specific significant features, so that deep features and shallow features are fused to achieve the purpose of feature fusion. Figure 4 and Figure 5 Shown

[0107] Furthermore, the loss function L is calculated as follows:

[0108] The loss function of the PD network includes the similarity loss L sim and smoothing loss L smooth Two parts; L sim The similarity loss is calculated by fixing the FA image and the deformed FA image, L smooth The smoothing loss is calculated by the deformation field; the calculation formula of the PD network loss function L is as follows:

[0109]

[0110] Among them, f and m are the three-dimensional space domain Ω∈R n In the image, f represents the fixed DTI, m represents the DTI to be registered, represents the deformation field, represents the deformed image output by the spatial transformation module, and ξ is the regularization parameter.

[0111] Since FA images clearly retain the organizational structure and position information of DTI, traditional registration methods based on DTI-derived images all use FA images to iteratively calculate the deformation field. Therefore, the main eigenvector is mainly used as complementary information to cooperate with FA images to obtain a more accurate deformation field, and does not participate in the calculation of the loss function.

[0112] Further, in step 3: the training method of the PD network is as follows:

[0113] The FA image of the DTI to be registered and the FA image of the fixed DTI in the training set are stacked as the input of the upper sub-network; the main eigenvector of the DTI to be registered and the main eigenvector of the fixed DTI in the training set are stacked as the input of the lower sub-network; the two sub-networks output the FA predicted deformation field and the main eigenvector predicted deformation field respectively, and the FA predicted deformation field and the main eigenvector predicted deformation field are input into the convolution layer together, and after the LeakyReLU activation function, a fused deformation field with 3 channels is obtained; the fused deformation field and the FA to be registered are output as the deformed FA after passing through the spatial transformation module.

[0114] Further, in step 4: the testing method of the PD network is as follows:

[0115] The FA images of the DTI to be registered and the fixed DTI in the test set are input into the upper sub-network of the trained PD network, and the main feature vectors of the DTI to be registered and the fixed DTI in the training set are input into the lower sub-network of the trained PD network; two deformation fields are obtained through the two sub-networks of the PD network, and the two deformation fields pass through a convolution layer and an activation function to output a fused deformation field; the fused deformation field and the DTI to be registered are input into the spatial transformation module to obtain the deformed DTI, and the deformed DTI is tensor redirected using the ANTs toolkit to finally obtain the registered DTI.

[0116] Furthermore, in step 5, the evaluation indexes of the registered DTI include subjective visual evaluation, registration time comparison, region-based evaluation and fiber bundle-based evaluation.

[0117] Further, in step 5, the evaluation method of the evaluation index is as follows:

[0118] The DTI output from step 4 and the fixed DTI selected from the IXI dataset were subjected to subjective visual evaluation, registration time comparison, region-based evaluation, and fiber bundle-based evaluation.

[0119] The region-based evaluation method includes two parts: tensor-based evaluation measure and scalar-based evaluation measure. The tensor-based evaluation measure and the scalar-based evaluation measure are compared through the data average value.

[0120] The tensor-based evaluation measures include tensor Euclidean distance and tensor deviation Euclidean distance; the scalar-based evaluation measures include bidirectional consistency, eigenvalue-eigenvector pair overlap rate, cross-correlation, and the angle between the main eigenvector pairs;

[0121] (I) Tensor Euclidean Distances (ED)

[0122] The smaller the ED value, the better the registration performance and the more aligned the biological tissue structure in DTI. The ED calculation formula is as follows:

[0123]

[0124] Where D1 and D2 represent two diffusion tensors respectively. The diffusion tensor is a 3×3 symmetric positive definite matrix. Tr() represents the trace of the tensor. The calculation formula of Tr() is as shown in formula (3):

[0125] Tr(D)=D xx +D yy +D zz (3)

[0126] Among them, D xx , D yy and D zz They are the three elements on the diagonal of the diffusion tensor D;

[0127] (II) Euclidean Distances of the Deviatoric tensor (DD)

[0128] The smaller the DD value, the better the registration performance and the more aligned the biological tissue structure in DTI. The DD calculation formula is as follows:

[0129]

[0130] Among them, ||·|| c Indicates calculation of Euclidean distance;

[0131] (III) Dyadic Coherence (κ)

[0132] The registration results are evaluated using the bidirectional consistency of the tensor. The changes in the biological tissue structure are evaluated by measuring the changes in the direction of the DTI main eigenvector. The κ calculation formula is as follows:

[0133]

[0134]

[0135] Among them, β i It is called a dyadic tensor; e1, e2, and e3 are the three eigenvectors of the diffusion tensor D; e1 is the main eigenvector; e 1x ,e 1y ,e 1z are the three components of the principal eigenvector respectively; is the principal eigenvector of the i-th subject; the larger the κ value, the better the performance of the registration algorithm;

[0136] (IV) Overlap of Eigenvalue-Eigenvector Pairs (OVL)

[0137] The larger the OVL value, the better the performance of the registration algorithm. The calculation formula of OVL is as follows:

[0138]

[0139] in, and are the jth eigenvalue and eigenvector of the i-th voxel of the registered DTI, respectively; and are the j-th eigenvalue and eigenvector of the i-th voxel of the fixed DTI, respectively;

[0140] 5. Cross Correlations (CC)

[0141] The cross-correlation between voxels is calculated using the scalar derivative measure of DTI. The larger the cross-correlation value, the higher the similarity of the images. The calculation formula of the cross-correlation CC is as follows:

[0142]

[0143] Where v represents the voxel of all images; X1(v) and X2(v) are scalar images; the value range of CC is [0,1];

[0144] The correlation between FA and Tr of the registered DTI and the fixed DTI was calculated respectively; CC Tr Calculation method and CC FA Same; FA's correlation CC FA The calculation formula is as follows (10):

[0145]

[0146] (VI) Angle between the principal eigenvector pairs

[0147] The angle (α) of the main eigenvector is:

[0148]

[0149] Among them, e′1 represents the main eigenvector of the fixed DTI, e1 represents the main eigenvector of the registered DTI, and the smaller the α value, the better the performance of the registration algorithm.

[0150] (VII) Evaluation based on fiber bundles: The fiber bundles of the registered DTI image are obtained by mapping the fiber bundles of the fixed image, and then the characteristic average is calculated and the average characteristic curve along the fiber bundle is drawn; the larger the normative correlation of the FA curves before and after registration, the better, the more similar the mean is, the better, and the smaller the variance is, the better; the FA characteristic curves of the selected fiber bundles include the genu of the corpus callosum (Genu), the splenium of the corpus callosum (Splenium), the left-superior thalamic radiation (L-str), the right-superior thalamic radiation (R-str), the left corticospinal / corticobulbar tracts (L-cst), and the right corticospinal / corticobulbar tracts (R-cst);

[0151] Taking the IXI dataset as an example, the registration results obtained by the dual-channel parallel DTI registration method (DTI RBPD) proposed in the present invention are compared and analyzed through multiple evaluation measures such as visual evaluation, registration time comparison, region-based evaluation, and fiber bundle-based evaluation. It can be concluded that the DTI registration method of DTI RBPD is better.

[0152] The DTI output from step 4 and the fixed DTI selected from the IXI dataset were used for subjective visual evaluation, registration time comparison, region-based evaluation, and fiber bundle-based evaluation; and compared with SyN, DTI-TK, and the DTI registration method based on VoxelMorph (Diffusion Tensor Image Registration method Based on VoxelMorph, i.e., DTI RBV).

[0153] Taking the IXI dataset as an example, the DTI images preprocessed by the IXI dataset were used to conduct training, testing, registration experiments, evaluation and comparative analysis on DTI RBPD.

[0154] 1. Visual evaluation

[0155] DTI RBPD has the best color preservation and is much higher than the other three methods in terms of the degree of detail correspondence of anatomical structures. Figures 6 to 8 shown. Figure 6 is the two-dimensional tangent direction of the axial direction, Figure 7is the two-dimensional tangent direction in the coronal direction, Figure 8 is the two-dimensional tangent direction in the sagittal position. In these three figures, sub-image (a) is the fixed DTI, (b) is the DTI to be registered, (c) is the registration result of the SyN method, (d) is the registration result of the DTI-TK method, (e) is the registration result of the DTI RBV method, and (f) is the registration result of the DTI RBPD method. Since the fixed and DTI to be registered are not brain data of the same individual, it is impossible for the results to be exactly the same. It is only necessary that the relative position of the anatomical structure is consistent and the tensor direction is roughly consistent.

[0156] Visually, the original DTI to be registered and the fixed DTI are obviously different. The DTI to be registered is larger overall. Figure 6 The red box in the axial position shows that compared with other comparison methods, the genu of the corpus callosum and the left and right suprathalamic radiations in the SyN registration results are significantly higher than those in other comparison methods. Figure 6 There are obvious differences in the relative positions of the anatomical structures (similar to the inverted triangle part in the figure). Figure 7 As can be seen from the red box in the coronal position, the SyN and DTI-TK registration results are difficult to distinguish in this part, and the anatomical structure is almost invisible. Although the DTI RBV and DTI RBPD methods can clearly see the anatomical structure in the box, it is difficult to tell which one is closer to the actual situation with the naked eye. Figure 8 As can be seen from the red box in the sagittal position, the DTI RBPD method has the best color preservation compared to the other three methods. At the same time, there is a thin green anatomical structure above the corpus callosum (red anatomical structure), and the DTI RBPD method corresponds best to the fixed image in terms of the spatial position of this structure. In addition, it can be seen from the yellow box that the registration result of the DTI-TK method is difficult to distinguish, while the anatomical structure of the DTIRBPD method is the best. In summary, the DTI RBPD method is much higher than the other three methods in terms of the degree of detail correspondence of the anatomical structure.

[0157] 2. Comparative evaluation of registration time consumption:

[0158] From the perspective of registration speed, the designed method is significantly faster than the traditional registration method and is basically the same as the DTI RBV method.

[0159] As shown in Table 1, DTI RBPD significantly improves the registration speed compared with traditional methods.

[0160] Table 1 Comparison of registration algorithm time consumption

[0161]

[0162] 3. Region-based evaluation

[0163] According to the numerical results in Table 2, DTI RBPD performs well in multiple measures and its registration performance is better than that of the comparison methods.

[0164] Table 2 Average values ​​of region-based evaluation metrics on the test set

[0165]

[0166]

[0167] 1) Tensor-based evaluation metrics

[0168] ED evaluates the registration effect by measuring the Euclidean distance between tensors. The smaller the ED value, the higher the similarity between tensors, the more accurate the structural alignment of brain white matter, and the better the registration algorithm performance. Table 2 shows the registration performance: DTI RBPD>DTI RBV>DTI-TK>SyN.

[0169] DD evaluates the registration effect by measuring the deviation Euclidean distance between tensors. The smaller the DD value, the higher the similarity between tensors, the more accurate the structural alignment of brain white matter, and the better the registration algorithm performance. Table 2 shows the registration performance: DTI RBPD>DTI RBV>DTI-TK>SyN.

[0170] 2) Scalar-based evaluation measures

[0171] κ focuses on the direction of the main eigenvector, thereby measuring the changes in brain white matter structure. The larger the value of κ, the higher the degree of symmetry of the eigenvector and the better the performance of the registration algorithm. Table 2 shows the registration performance: DTI RBPD>DTI RBV>DTI-TK>SyN.

[0172] OVL measures the overlap rate of each pair of eigenvalues ​​and eigenvectors of the fixed image and the registered image tensor. The larger the OVL value, the higher the overlap rate and the better the registration algorithm performance. Table 2 shows the registration performance: DTI-TK>DTI RBV>DTI RBPD>SyN.

[0173] The present invention calculates the cross-correlation of FA value and Tr value between the fixed image and the registered image respectively. The larger the value, the higher the image similarity and the better the registration algorithm performance. Table 2 shows the registration performance: DTI RBPD>DTI RBV>DTI-TK>SyN.

[0174] The degree of alignment between the two images can be measured by measuring the angle α between the main eigenvectors of the fixed DTI and the registered DTI. The smaller the angle, the higher the degree of alignment and the better the registration algorithm performance. Table 2 shows the registration performance: DTI-TK>DTI RBV>DTI RBPD>SyN.

[0175] 4. Evaluation based on fiber bundles

[0176] According to the canonical correlation of FA characteristics between the fiber bundles of the registered DTI and the fixed DTI (as shown in Table 3), and the number of registration failures of different methods in different regions of interest (ROI) (as shown in Table 4), it is concluded that the registration performance of DTI RBPD is better than that of the comparison methods;

[0177] The FA characteristic curves of the selected fiber bundles include the Genu of the corpus callosum (Genu), the Splenium of the corpus callosum (Splenium), the left-superior thalamic radiation (L-str), the right-superior thalamic radiation (R-str), the left corticospinal / corticobulbar tracts (L-cst), and the right corticospinal / corticobulbar tracts (R-cst).

[0178] Table 3 Canonical correlation of FA properties between the registered DTI fiber bundles and the fixed DTI fiber bundles

[0179]

[0180] Note: * indicates significant correlation at the 0.01 level

[0181] Table 3 shows the registration performance: DTI RBPD>DTI RBV>DTI-TK>SyN. Although DTI RBPD in Table 3 performs poorly on the two ROIs of Splenium and L-str, the correlation between the proposed method and the optimal algorithm DTI-TK is relatively small. The evaluation based on the fiber bundle is performed on a specific ROI, not on the DTI individual in the entire space. Therefore, poor performance on individual ROIs does not affect the evaluation of the overall performance. In general, the registration performance of DTI RBPD is better than that of the comparison method.

[0182] Table 4 Number of registration failures of four registration methods on different ROIs

[0183]

[0184] Based on the correlation coefficient, the present invention also calculates the number of subject ROIs that are not correctly registered to the fixed image. The threshold is selected as 0.85, and ROIs with a correlation coefficient lower than 0.85 are marked as registration failures. Four methods are used to register the test individuals (20), and the number of registration failures on different ROIs is shown in Table 4. It can be seen from the table that DTI RBPD has the least number of failures and the best registration performance.

Claims

1. A DTI registration method based on a dual-channel parallel network, characterized in that: The DTI registration process is as follows: Step 1, determine the data set: select the IXI data set, preprocess the IXI data set to obtain the DTI data set, and divide the DTI data set into a training set and a test set; Step 2: Build a PD network; Step 3, training the PD network; Step 4, test the PD network; Step 5, evaluate the registered DTI; Step 6: Verify the generalizability of the DTI registration method on different datasets: Use the dHCP dataset and repeat steps 1 to 5 for verification.

2. A DTI registration method based on a dual-channel parallel network according to claim 1, characterized in that: In step 1, the IXI dataset preprocessing method is as follows: First, the diffusion image in the IXI dataset is processed by eddy current correction to obtain Image 1; then the skull in Image 1 is removed to obtain Image 2; finally, a DTI dataset with an individual resolution of 128×128×64 is obtained by tensor calculation; the DTI dataset is divided into a training set and a test set in a ratio of 8:1; A DTI with excellent morphology and image quality was manually observed and randomly selected from the DTI dataset as the fixed DTI in the registration process, and the remaining DTIs were the DTIs to be registered. The derived FA images and main eigenvectors of each DTI in the DTI dataset were calculated using FSL software. The FA image included the fixed FA and the FA to be registered. The main eigenvector included the fixed main eigenvector and the main eigenvector to be registered. The FA image pair was the fixed FA and the FA to be registered, and the main eigenvector pair was the fixed main eigenvector and the main eigenvector to be registered.

3. A DTI registration method based on a dual-channel parallel network according to claim 1, characterized in that: In step 2, the PD network construction method is as follows: The PD network includes two identical sub-networks and a spatial transformation module; the sub-network is a VoxelMorph network that introduces an attention mechanism; the sub-network for inputting FA image pairs is the upper sub-network, and the sub-network for inputting main feature vector pairs is the lower sub-network; S2.1, input the FA image pair into the upper sub-network, and input the main feature vector pair of the DTI to be registered and the fixed DTI into the lower sub-network; the FA image and DTI are output as two deformation fields through the two sub-networks respectively; S2.2, the two deformation fields pass through a convolution layer and an activation function, and the output is a fused deformation field, and the number of output channels of the convolution layer is 3; S2.3, fusion of the deformation field and the DTI to be registered into the spatial transformation module to obtain the deformed DTI; S2.4, the deformed DTI is output as the registered DTI after tensor redirection; S2.5, calculate the loss function L of the PD network.

4. A DTI registration method based on a dual-channel parallel network according to claim 3, characterized in that: The sub-network includes an encoder, a series skip connection and a decoder; The encoder uses a three-dimensional convolutional layer with a convolution kernel size of 3×3×3 and a step size of 2, and the convolutional layer is followed by a LeakyReLU activation function; the encoder reduces the spatial dimension of the feature map through convolution to obtain DTI shallow features; The decoder uses a three-dimensional convolutional layer with a convolution kernel size of 3×3×3 and a step size of 2. The convolutional layer is followed by a LeakyReLU activation function. The decoder enlarges the spatial dimension of the feature map through upsampling and convolutional layers to obtain DTI deep features. The spatial dimension of the decoder changes in the opposite direction to that of the encoder. The number of output channels of the decoder is 3, which corresponds to the number of channels of the predicted deformation field. The serial jump connection is a CBAM that connects channel attention and spatial attention in series. CBAM assigns weights to deep features and shallow features to fuse the deep features with the shallow features.

5. The DTI registration method based on a dual-channel parallel network according to claim 3, characterized in that: The calculation method of the loss function L is as follows: The loss function of the PD network includes the similarity loss L sim and smoothing loss L smooth Two parts; L sim The similarity loss is calculated by fixing the FA image and the deformed FA image, L smooth The smoothing loss is calculated by the deformation field; the calculation formula of the PD network loss function L is as follows: Among them, f and m are the three-dimensional space domain Ω∈R n In the image, f represents the fixed DTI, m represents the DTI to be registered, represents the deformation field, represents the deformed image output by the spatial transformation module, and ξ is the regularization parameter.

6. The DTI registration method based on a dual-channel parallel network according to claim 1, characterized in that: In step 3: the training method of the PD network is as follows: The FA image and main eigenvector in the training set are input into the PD network; the two deformation fields of the PD network are input into the convolution layer together, and after the LeakyReLU activation function, a fused deformation field with 3 channels is obtained; the fused deformation field and the FA to be registered are output as the deformed FA after passing through the spatial transformation module.

7. The DTI registration method based on a dual-channel parallel network according to claim 1, characterized in that: In step 4: the test method of the PD network is as follows: The FA images in the test set and the input main eigenvector are input into the trained PD network; two deformation fields are obtained through the PD network, and the two deformation fields pass through a convolution layer and an activation function at the same time to output a fused deformation field; the fused deformation field and the DTI to be registered are input into the spatial transformation module to obtain the deformed DTI; the deformed DTI uses the ANTs toolkit for tensor redirection, and finally the registered DTI is obtained.

8. The DTI registration method based on a dual-channel parallel network according to claim 1, characterized in that: In step 5, the evaluation indexes of the registered DTI include subjective visual evaluation, registration time comparison, region-based evaluation and fiber bundle-based evaluation.

9. The DTI registration method based on a dual-channel parallel network according to claim 1, characterized in that: In step 5, the evaluation method of the evaluation index is as follows: The region-based evaluation method includes two parts: tensor-based evaluation measure and scalar-based evaluation measure. The tensor-based evaluation measure and the scalar-based evaluation measure are compared through the data average value. The tensor-based evaluation measures include tensor Euclidean distance and tensor deviation Euclidean distance; the scalar-based evaluation measures include bidirectional consistency, eigenvalue-eigenvector pair overlap rate, cross-correlation, and the angle between the main eigenvector pairs; The calculation formula of tensor Euclidean distance ED is as shown in formula (2): Where D1 and D2 represent two diffusion tensors respectively. The diffusion tensor is a 3×3 symmetric positive definite matrix. Tr() represents the trace of the tensor. The calculation formula of Tr() is as shown in formula (3): Tr(D)=D xx +D yy +D zz (3) Among them, D xx , D yy and D zz They are the three elements on the diagonal of the diffusion tensor D; The calculation formula of tensor deviation Euclidean distance DD is as shown in formula (4): Among them, ||·|| c Indicates calculation of Euclidean distance; The calculation formula of bidirectional consistency κ is as follows: Among them, β i It is called a dyadic tensor; e1, e2, and e3 are the three eigenvectors of the diffusion tensor D; e1 is the main eigenvector; e 1x ,e 1y ,e 1z are the three components of the principal eigenvector respectively; is the principal eigenvector of the i-th subject; the larger the κ value, the better the performance of the registration algorithm; The calculation formula of the eigenvalue-eigenvector overlap ratio OVL is as shown in formula (8): in, and are the jth eigenvalue and eigenvector of the i-th voxel of the registered DTI, respectively; and are the j-th eigenvalue and eigenvector of the i-th voxel of the fixed DTI, respectively; The calculation formula of the cross-correlation CC is as follows: Where v represents the voxel of all images; X1(v) and X2(v) are scalar images; the value range of CC is [0,1]; The correlation between FA and Tr of the registered DTI and the fixed DTI was calculated respectively; CC Tr Calculation method and CC FA Same; FA's correlation CC FA The calculation formula is as follows (10): The angle (α) between the principal eigenvector pairs is: Among them, e′1 represents the main eigenvector of the fixed DTI, e1 represents the main eigenvector of the registered DTI, and the smaller the α value, the better the performance of the registration algorithm; Evaluation based on fiber bundles: The fiber bundles of the registered DTI image are obtained by mapping the fiber bundles of the fixed image, and then the average characteristic value is calculated and the average characteristic curve along the fiber bundle is plotted.

Images