Multi-modal medical image intracranial aneurysm segmentation method based on geometric deep learning

Abstract

The invention relates to a multi-modal medical image intracranial aneurysm segmentation method based on geometric deep learning, and the method comprises the steps: 1, designing a medical diffusion network based on a deep neural network of frequency domain convolution; 2, preprocessing the three-dimensional grid data to realize the manifold of the model and calculate the characteristics of the grid model; and step 3, using a multiple loss function to reduce the influence of the unbalanced data set, and the like. According to the method, the tedious process of manual feature design and multi-stage processing in a traditional segmentation method is avoided; the MedDiffusion Net designed by the method solves the problem of intracranial aneurysm segmentation of medical images of various modalities, and experimental results on MRA, CTA and DSA data show that the method can obtain accurate intracranial aneurysm segmentation without training data of different modalities respectively; and adaptively fusing the features in the neighborhood by learning a diffusion parameter t.

Description

Technical Field

[0001] The present invention belongs to the technical field of geometric deep learning, and in particular relates to a method for segmenting intracranial aneurysms in multimodal medical images based on geometric deep learning. Background Art

[0002] Intracranial aneurysm is a relatively common high-risk cerebrovascular disease involving localized pathological dilatation of intracranial arteries. Over time, the development and expansion of aneurysms increase the risk of subarachnoid hemorrhage (SAH) and hemorrhagic stroke. The rupture of intracranial aneurysms and their direct complications lead to the death of more than 25% of patients, and more than 50% of patients suffer from permanent neurological damage. The prevalence of intracranial aneurysms in the natural population is 3% to 8%, and the overall annual risk of rupture is 0.95% to 2%, while the risks of endovascular treatment and craniotomy and clipping are 5.3% and 6.3%, respectively.

[0003] In clinical practice, the diagnosis of intracranial aneurysms relies on the analysis of three-dimensional medical imaging data, including computed tomography angiography (CTA), three-dimensional rotational digital subtraction angiography (3D-DSA), and magnetic resonance angiography (MRA). At present, doctors mainly manually annotate aneurysms through the two-dimensional views of these 3D brain images, which is labor-intensive, time-consuming, and easy to miss lesions. In addition, due to the hemodynamic effects and the complexity and variability of cerebrovascular structures, aneurysms show significant individual differences in morphology, size, and number, which further increases the difficulty of segmentation. In order to improve segmentation efficiency and accuracy and reduce the potential risk of aneurysm rupture, it is necessary to develop efficient and accurate aneurysm segmentation methods to greatly reduce the risk of aneurysm rupture.

[0004] With the significant research in the field of computer vision, deep learning has been widely used in medical image analysis. Convolutional neural network (CNN) is a deep learning architecture. CNN consists of multiple convolutional layers and pooling layers, which can extract multi-level features of the image layer by layer. Shallow convolutional layers extract low-level features, such as the edge of an aneurysm, and deep convolutional layers can extract high-level features, such as the shape of an aneurysm. This multi-level feature extraction capability enables CNN to better understand and segment complex medical images.

[0005] U-Net is a typical CNN-based medical image segmentation network with a U-shaped structure, including an encoder and a decoder. The encoder extracts features layer by layer, the decoder restores the spatial resolution layer by layer, and passes the features in the encoder directly to the corresponding decoder layer through skip connections to retain more detailed information. U-Net performs well in various medical image segmentation tasks and is widely used in tumor segmentation, organ segmentation and other fields. However, most image-based aneurysm segmentation methods are only applicable to one modality of medical images (CTA, DSA, TOF-MRA, etc.), which limits the practical application of such methods. Summary of the invention

[0006] Based on the above defects of the prior art, the present invention proposes a multimodal medical image intracranial aneurysm segmentation method based on geometric deep learning.

[0007] The method for segmenting intracranial aneurysms in multimodal medical images based on geometric deep learning of the present invention comprises the following specific steps:

[0008] Step 1: Design the Medical Diffusion Network (MedDiffusionNet) based on the deep neural network of frequency domain convolution

[0009] Step 1.1, Instance Normalization

[0010] Instance normalization (IN) normalizes the features of each channel to satisfy the normal distribution. IN normalizes each sample individually, while BN normalizes all samples in the batch. The specific calculation method of IN is as follows:

[0011]

[0012] Among them, x ic represents the feature on channel c of the i-th sample without instance normalization, then y ic represents the feature on channel c of the ith sample after instance normalization. γ and β are learnable parameters in instance normalization. Their function is to linearly transform the features that have been normalized to the standard normal distribution so as not to limit the expression ability of the network. μ(x ic ) represents x ic The expectation of σ(x ic ) represents x ic The standard deviation of

[0013] Step 1.2, establish a feature linear transformation layer;

[0014] The feature linear transformation layer is used to achieve the change of feature dimension through a learnable linear transformation. The feature linear transformation layer is designed in the front-end module of the network to increase the feature dimension and match the hidden layer scale, while the feature linear transformation layer in the back-end of the network is to reduce the feature dimension and match the classification category. The specific calculation of the feature linear transformation layer is as follows (2):

[0015] x′=xW ......(2),

[0016] Where x represents a sample input feature, x′ represents the feature of the sample after feature linear transformation, W represents the matrix for feature linear transformation, and W is also the only parameter to be learned in the feature linear transformation layer. The calculation method of the feature linear transformation of the following medical diffusion network is consistent with the above formula (1), and the difference lies in the dimension, especially the dimension of W;

[0017] Step 1.3, frequency domain convolution based on heat diffusion

[0018] The heat diffusion equation describes the diffusion of heat in a certain space. When this equation is extended to the non-Euclidean domain, the specific equation is as follows (3):

[0019]

[0020] Where Δ represents the Laplace–Beltrami operator (Laplace operator for short); u(t) represents the function of the heat at each location in the model space at time t. In the case of discreteness, the LB operator Δ is simplified to the form of an LB matrix, and u(t) is simplified from a function to a finite-dimensional vector. The heat diffusion equation (3) can be regarded as a system of first-order linear homogeneous differential equations with constant coefficients. When the boundary condition is t = 0 and u(0) is known, the solution of u(t) is as follows (4):

[0021] u(t)=e -tΔ u(0) ......(4),

[0022] Among them, e -tΔ The matrix index of e is as follows:

[0023]

[0024] After performing eigenvalue decomposition on the matrix A, equation (5) is further simplified. The LB matrix Δ is a real symmetric matrix. Perform eigenvalue decomposition on it. Let the matrix composed of the eigenvectors of the Laplace matrix be U, and the diagonal matrix composed of its corresponding eigenvalues ​​be Λ. Then the Laplace matrix is ​​transformed into Δ = UΛU T , as above, we get the following formula (6):

[0025] u(t)=Ue -tΛU T u(0) ......(6),

[0026] Where -tΛ is a diagonal matrix, and the matrix exponentiation of the diagonal matrix e is converted to the direct element exponentiation to simplify the u(t) solution process;

[0027] Convolution is defined by the heat diffusion equation as shown above. Convolution can realize feature fusion in the spatial domain and feature fusion between channels. The heat diffusion equation (6) above is used to realize feature fusion in the spatial domain. The initial feature is set to u(0). The diffusion time t is the parameter to be learned. The size of the fusion area is implicitly determined by t. When t is small, the features in a small area are fused; when t is large, the features in a larger range are fused; when t tends to infinity, the features of all areas in the space are fused. This is similar to the definition of frequency domain convolution. -tΛ The convolution kernel is regarded as a frequency domain, and the whole process is regarded as a convolution in the frequency domain;

[0028] Step 1.4: Establish a learnable gradient feature module

[0029] The heat diffusion equation is isotropic, which limits the expressive power of the convolution operator. The gradient of the calculated features is introduced to increase its expressive power and realize anisotropic feature fusion, as follows:

[0030] In order to calculate the gradient of the function defined on the mesh vertex domain, the function on the mesh vertex domain is regarded as a simplified representation of the function defined on the mesh manifold (only the function value at each vertex is given, and the function values ​​at other positions are obtained by interpolation of the function values ​​of the surrounding vertices). For the function defined on the manifold, the gradient at any point is a two-dimensional vector located in its tangent plane, and the two-dimensional vector points to the direction where the function value changes the most. For convenience, complex numbers are used to represent the gradient at any point. For the function on the mesh manifold, the gradient at each point on a face is the same, and the gradient between different faces is different. The gradient on the edge and vertex is approximated by the following method:

[0031] First, the neighborhood vertices are projected onto the current point tangent plane;

[0032] Then the function is Taylor expanded at the current point to obtain a set of gradient-related equations;

[0033] Finally, the least squares method is used to solve the equations. The gradient on the vertex is quickly obtained using the pre-calculated gradient matrix G. After the gradient function of each feature function is calculated, the learnable gradient feature is calculated according to the following formula (7):

[0034] h = tanh(<z,Az> ) ......(7), where z i ∈C HSrepresents the gradient corresponding to all hidden layer features at grid vertex i; A represents a parameter matrix to be learned. The parameter matrix can be a real matrix or a complex matrix. Complex matrices are used in subsequent experiments; <·, ·> represents the calculation of the dot product of the complex numbers at corresponding positions in two complex vectors. The dot product here refers to the dot product of the two-dimensional vector represented by the complex numbers; tanh is an activation function that can ensure that the final feature range is [-1,1]; the final calculation result h i ∈R HS It is the learnable gradient feature at vertex i;

[0035] The learnable gradient feature is to first perform a learnable linear combination on the gradient of the input feature at each grid vertex, then calculate the dot product between the result of the gradient linear combination and the original gradient, and finally use the activation function to process the result of the dot product and use it as the final learnable gradient feature;

[0036] Since the solution of the heat diffusion equation is calculated in the frequency domain, in order to simplify and unify the calculation process, it is necessary to calculate the gradient function of the characteristic function in the frequency domain. To this end, it is only necessary to do some processing on the gradient matrix G, let Where U represents the matrix formed by concatenating the eigenvectors of the LB matrix. Then we only need to By multiplying it with the function transformed to the frequency domain, we can get the gradient function of the function in the spatial domain;

[0037] Step 1.5, Multilayer Perceptron Module

[0038] In order to achieve the fusion of feature information between different channels, a multi-layer perceptron (MLP) is added at the end of each diffusion block. The input of MLP consists of three parts, namely, the original input feature x, the gradient feature obtained by the learnable gradient feature module, and the diffusion result returned to the spatial domain by the diffusion vector calculated in the frequency domain;

[0039] Step 2: Preprocess the 3D mesh data to realize the manifold of the model and calculate the features of the mesh model. It is necessary to sequentially process the 3D mesh data to realize the manifold of the mesh model, calculate the features of the mesh model, and calculate the Laplace-Beltrami matrix:

[0040] Step 2.1, processing the reconstructed mesh model to make it a manifold;

[0041] Step 2.2, calculate the point and surface features of the three-dimensional mesh model. The Euclidean space coordinates, wave kernel signature (WKS), maximum principal curvature, minimum principal curvature, average geodesic distance (AGD), shape diameter function (SDF) and other features are calculated on the vertices of the mesh model for subsequent learning, where:

[0042] Euclidean space coordinates are the basic features of the grid model, but in order to train together with other features, they need to be normalized so that each dimension of their three-dimensional coordinates is in the range of [-1,1];

[0043] The wave nuclear signature WKS comes from the theory of spectral analysis and is an advanced shape feature. WKS is derived from the Schrödinger equation and describes the average probability distribution of quantum particles of different energy levels at each vertex on the grid. WKS will not change due to rigid body transformations such as translation and rotation of the grid and is a typical intrinsic feature.

[0044] The maximum principal curvature refers to the maximum curvature value among all possible curves passing through a point, which describes the direction in which the surface is most sharply curved at that point. The minimum principal curvature refers to the minimum curvature value among all possible curves passing through that point, which describes the direction in which the surface is most gently curved at that point.

[0045] The average geodesic distance AGD describes the average geodesic distance from a vertex on the grid to all other points on the grid. When AGD is larger, it means that the point is more likely to be at the end of the grid model; when AGD is smaller, the point is located at the center of the grid model.

[0046] The shape diameter function SDF defines a scalar value at each point of the three-dimensional grid. This value is obtained by calculating the distance to the opposite side of the shape along the internal direction of the shape (i.e., the normal vector of the surface) at the point and averaging it in the neighborhood of the point. The value of the shape diameter function geometrically represents the thickness of the local area of ​​the shape, thereby describing the overall structure and local characteristics of the blood vessel or aneurysm.

[0047] Step 2.3, calculate the Laplace-Beltrami (LB) matrix on the triangular mesh model. The Laplace operator is a differential operator for functions whose domain is Euclidean space. After being acted on, the function will become the divergence of its own gradient. The Laplace-Beltrami operator is also a differential operator for functions whose domain is a manifold. Similar to the Laplace operator, the function after being acted on will also become the divergence of its own gradient. When facing a mesh manifold, the continuous function on it can be simplified by the multidimensional vector of the function value at each vertex of the mesh (the function value at other positions can be obtained by interpolation of the function value of the surrounding vertices). At this time, the LB operator for the function defined on the domain also degenerates into the Laplace-Beltrami matrix (LB matrix). The LB matrix is ​​calculated by using the cotan algorithm. The result obtained by the algorithm is called the cotan-Laplace matrix. The calculation method of the cotan-Laplace matrix when acting on any function is as follows (8):

[0048]

[0049] Where Δ represents the Cotan-Laplace matrix, u represents a function defined on the grid manifold, (Δu) i It represents the value of u at vertex i after Δ acts on it, that is, the divergence of the gradient of u at vertex i; represents the set of all neighboring vertices of vertex i; α j and β j represents the diagonal angle between vertex i and adjacent vertex j, u i and u j represents the function value of function u at vertex i and vertex j. The above formula describes the calculation method of the Cotan-Laplacian matrix when it acts on a function defined on a mesh manifold. From it, the specific form of the Cotan-Laplacian matrix can be derived as follows:

[0050]

[0051] Among them, Δ i，jrepresents the value at the i-th row and j-th column of the Cotan-Laplacian matrix; the other symbols are consistent with those in formula (8). It can be seen that the Cotan-Laplacian matrix is ​​a symmetric matrix. According to the properties of the symmetric matrix, the Cotan-Laplacian matrix can be decomposed into eigenvalues ​​to obtain its eigenvalues ​​and the corresponding eigenvectors. After obtaining the eigenvalues ​​and eigenvectors of the LB matrix, the frequency domain can be defined on the grid manifold in the same way as the frequency domain in the Euclidean space. For a function whose domain is the Euclidean space, the airspace refers to the Euclidean space where the domain is located, and the frequency domain refers to a functional space formed by the eigenfunction of the Laplace operator as the basis. When it is extended to the function on the grid manifold, the airspace refers to the grid manifold, and the frequency domain refers to the space formed by the eigenvectors of the LB matrix.

[0052] Step 3: Use multiple loss functions to reduce the impact of unbalanced datasets

[0053] The multiple loss function consists of three parts, namely, BCE loss function, Focal loss function and Dice loss function. The multiple loss function effectively combines the advantages of the above three loss functions and improves the accuracy of intracranial aneurysm segmentation, as shown in the following formula (10):

[0054]

[0055] The BCE loss function effectively quantifies the difference between the predicted segmentation results and the actual segmentation results, providing a robust metric for the classification tasks of blood vessels and aneurysms.

[0056]

[0057] in, and l i Represent the predicted label and true label of the i-th sample respectively, and N represents the batch size during training. In addition, the Focal loss function enhances the segmentation performance of intracranial aneurysms by adjusting the weights of simple samples and difficult samples, as well as the weights of the vascular part and the aneurysm part, as shown in the following formula (12):

[0058]

[0059] in, l i and N are defined the same as in the BCE loss function. The weighting factor α balances the weights of vessels / aneurysms, while the focusing factor γ controls the balance of easy / difficult examples. α and γ are specified as 0.8 and 1.5, respectively.

[0060] Finally, the Dice loss function alleviates the impact of the imbalanced ratio of blood vessels and aneurysms, while enhancing the sensitivity at their segmentation boundaries. The Dice loss function is shown in the following formula (13):

[0061]

[0062] in, and l i Same definition as formula (12).

[0063] The method of the present invention has the following beneficial effects:

[0064] 1. The method of the present invention designs MedDiffusionNet by improving the frequency domain neural network. The present invention avoids the tedious process of manually designing features and multi-stage processing required by traditional segmentation methods.

[0065] 2. MedDiffusionNet designed by the method described in the present invention solves the problem of intracranial aneurysm segmentation in multi-modal medical images. Experimental results on MRA, CTA and DSA data show that the method described in the present invention can obtain accurate intracranial aneurysm segmentation results without separate training of data of different modalities.

[0066] 3. The method of the present invention realizes feature fusion in the spatial domain and feature fusion between channels based on the frequency domain convolution of the heat diffusion equation, and adaptively fuses features within the neighborhood by learning the diffusion parameter t.

[0067] 4. The MedDiffusionNet designed by the method described in the present invention adds instance normalization to improve the convolution feature extraction effect based on the heat diffusion equation in the network, reduce the risk of overfitting, and ensure that the network obtains better results.

[0068] 5. The MedDiffusionNet designed by the method described in the present invention adopts multiple loss functions, which can adjust the training weights of the blood vessel part and the aneurysm part and the importance of simple samples and difficult samples, thereby alleviating the differences in the mesh models reconstructed from images of different modalities.

[0069] 6. The method of the present invention realizes the automatic quantitative calculation and analysis of intracranial aneurysms in 3D grid morphology, including aneurysm height, aneurysm neck width, aneurysm volume and other morphological indicators. BRIEF DESCRIPTION OF THE DRAWINGS

[0070] Figure 1 (a) to (c) show schematic diagrams of the MedDiffusionNet network structure of the method of the present invention;

[0071] Figures 2 to 4It is a schematic diagram of the intracranial aneurysm segmentation result of the method of the present invention on a mixed data set;

[0072] Figure 5 (a) to (e) are schematic diagrams of aneurysm morphology indicators used in the method of the present invention;

[0073] Figure 6 This is a schematic diagram of the error estimation between the intracranial aneurysm segmentation result and the actual aneurysm;

[0074] Figure 7 (a) to (f) are schematic diagrams of consistency analysis between the intracranial aneurysm segmentation results of the present invention and the real aneurysm. DETAILED DESCRIPTION

[0075] The following is combined with the instructions attached Figures 1 to 7 The technical solution of the method of the present invention is further supplemented by the following embodiments.

[0076] Example

[0077] like Figure 1 As shown, the embodiment of the method of the present invention mainly includes:

[0078] Step 1 designs a medical diffusion network based on a deep neural network of frequency domain convolution, wherein the collection and processing of data sets uses a mixed data set composed of three data sets as a verification of the invention content, including the IntrA data set, the IntrANeurIST data set and the IntrANeurCT data set, which are reconstructed from the original medical images of MRA, DSA and CTA respectively, and each sample is a blood vessel-aneurysm mesh model, wherein the IntrA data set and the IntrANeurIST data set are public data sets, which contain the annotation information of blood vessels and aneurysms, the IntrA data set contains 116 samples, and the IntrANeurIST data set contains 56 samples. Compared with the samples of the IntrA data set, the proportion of aneurysms in the blood vessel-aneurysm model of the data set is relatively small, which also increases the difficulty of learning the network model; IntrANeurCT contains 8 samples, which are derived from the CTA images of five cases. The original CTA data are subjected to threshold segmentation processing by the VMTK tool, the vascular body data therein is extracted and converted into a three-dimensional mesh model, and finally the aneurysm part is manually annotated;

[0079] Since the mesh models in the IntrANeurIST dataset contain several times the number of meshes as the IntrA dataset, and the triangular meshes in the IntrANeurIST dataset are highly irregular, in contrast to the relatively regular shapes of the triangles in the IntrA dataset, this leads to an imbalance in the samples within the dataset. In addition, the two ends of the blood vessels in the IntrA dataset are non-closed, while the samples in the IntrANeurIST dataset and the IntrANeurCT dataset are all closed models. In order to mitigate the impact of these factors during training, the method of the present invention closes the holes in 116 samples of the IntrA dataset in the mixed dataset, ensuring that each mesh model has a closed surface after repair, and performs a re-meshing operation on each of the 180 mesh models, standardizing each blood vessel-aneurysm model to approximately 4,000 vertices.

[0080] Step 2 preprocesses the three-dimensional mesh data, wherein the calculation of mesh model related features is divided into two parts, firstly, the calculation of vertex geometric features, for all experimental data, MATLAB is used to calculate the 21-dimensional WKS features, maximum principal curvature, minimum principal curvature and AGD on each vertex, CGAL is used to calculate the SDF on each face, and then converted into the SDF on each vertex, in the method described in the present invention, 6-dimensional WKS features, maximum principal curvature, minimum principal curvature, AGD, SDF and vertex space features such as 3-dimensional Euclidean space coordinates are selected, totaling 13-dimensional features, and then Numpy and Scipy are used to calculate the LB matrix Δ, gradient matrix G and LB matrix eigenvalues ​​Λ and corresponding eigenvectors U of each mesh model. All the above features are stored on the hard disk to avoid repeated calculation of these items during the training process, thereby reducing the network training time, and the method described in the present invention divides the data set and performs subsequent training through a five-fold cross-validation strategy, each fold contains 150 training set samples and 30 test set samples, and the test set samples of each fold will not appear repeatedly. The vertex geometric features calculated in (2) are used as the input of the network model. When implementing the method described in the present invention, the network contains 8 stacked diffusion blocks, the size of the hidden layer in all diffusion blocks is set to 512, the initial learning rate is set to 0.0005, and the Adam optimizer and cosine annealing strategy are used to adjust the learning rate. Each fold is trained for 300 rounds.

[0081] Step 3 uses multiple loss functions to reduce the impact of unbalanced data sets. The multiple loss functions involved in the method of the present invention are composed of three parts, namely BCE loss function, Focal loss function and Dice loss function. The multiple loss function effectively combines the advantages of the above three loss functions, improves the accuracy of intracranial aneurysm segmentation, and uses the trained neural network model to segment the blood vessel-aneurysm model. First, the relevant features of the model to be segmented are calculated, and then they are input into the trained MedDiffusionNet model. Finally, the output prediction label is obtained to complete the segmentation of the three-dimensional mesh model. Figures 2 to 4 The segmentation result of the test set for sampling five-fold cross validation;

[0082] The method of the present invention is evaluated by calculating the morphological indexes of the aneurysm segmentation results. Figure 5 The definitions of some morphological indicators are shown, and the specific calculation methods are as follows: Neck width (NW) is defined as the distance between the two farthest points of the aneurysm neck (for example, point A and point B); Neck area (NA) is the area defined by the vertices on the aneurysm neck projected on the plane, and the area of ​​the projection plane is calculated by summing the areas of many triangles, each triangle is delimited by the projection center point and the boundary points; Aneurysm volume (AV) is calculated by summing the volumes of many tetrahedrons, each tetrahedron is defined by the center point of the aneurysm neck and all the faces on the aneurysm; Aneurysm area (AA) is defined as the outer surface area excluding the neck area; Aneurysm depth (AD) is defined as the maximum distance from the aneurysm vertex to the center of the aneurysm neck (for example, point C and point D); Aspect ratio (AR) is defined as the ratio of aneurysm depth to neck width;

[0083] The above morphological indicators can reflect the morphology of aneurysms. By comparing with the actual aneurysm morphological indicators, the effectiveness of the method of the present invention can be verified. Figure 6 and Figure 7 It is the result of morphological analysis of predicted aneurysms and real aneurysms. Figure 6 The violin plot used shows the error estimation between the predicted aneurysm and the real aneurysm. The method of the present invention uses a threshold of 1.5 times the interquartile range (IQR) to identify outliers in the data. There are three dotted lines in each violin plot, the upper line represents the third quartile, the middle line represents the second quartile, and the lower line represents the first quartile. It can be observed that the difference between the predicted segmentation result of the method of the present invention and the real value is small. Since the method of the present invention uses the least squares method to fit the projection plane and calculate its area, the change of the aneurysm neck vertex will cause significant fluctuations, which further leads to a relatively large error in the measurement of the neck area. Figure 7The Altman graph is used to show the morphological consistency analysis between the intracranial aneurysm segmentation results and the real aneurysm. Each point in the graph corresponds to the average and difference of the measurement pair (predicted aneurysm and real aneurysm), the average of the two measurements is used as the x-coordinate, and the difference between these values ​​is used as the y-coordinate. If most points are within the consistency range (depending on the morphological indicators of the real aneurysm), it indicates that there is a strong consistency between the two sets of results; in addition, if the average deviation is zero, it indicates that there is no systematic deviation between them. According to the drawn graph, the segmentation results of the method described in the present invention are highly consistent with the real situation, the average difference of all metrics is close to zero, and the vast majority of samples are within the consistency range.

[0084] The above description is only a specific implementation mode of the present invention, but the protection scope of the present invention is not limited thereto. Any changes or substitutions that can be easily thought of by any technician familiar with the technical field within the scope disclosed in the present invention should be covered by the protection scope of the claims of the present invention.

Claims

1. A method for segmenting intracranial aneurysms in multimodal medical images based on geometric deep learning, comprising the following steps: Step 1, design a medical diffusion network based on a deep neural network with frequency domain convolution; Step 2, preprocessing the three-dimensional grid data, realizing the manifold of the model and calculating the characteristics of the grid model, it is necessary to sequentially perform the manifold processing of the grid model, calculate the characteristics of the grid model and calculate the Laplace-Beltrami matrix on the three-dimensional grid data; Step 3: Use multiple loss functions to reduce the impact of imbalanced datasets.

2. According to the multimodal medical image intracranial aneurysm segmentation method based on geometric deep learning as claimed in claim 1, the deep neural network based on frequency domain convolution in step 1 is used to design a medical diffusion network: Step 1.1, Instance Normalization The features of each channel are normalized to satisfy the normal distribution. IN is the normalization of each sample individually, while BN is the uniform normalization of all samples in the batch. The specific calculation method of IN is as follows: in, represents the feature on channel c of the i-th sample without instance normalization, then y ic represents the feature on channel c of the ith sample after instance normalization. γ and β are learnable parameters in instance normalization. Their function is to linearly transform the features that have been normalized to the standard normal distribution so as not to limit the network's expressive power. μ(x ic ) represents x ic The expectation of σ(x ic ) represents x ic The standard deviation of Step 1.2, establish a feature linear transformation layer; The feature linear transformation layer is used to achieve the change of feature dimension through a learnable linear transformation. The feature linear transformation layer is designed in the front-end module of the network to increase the feature dimension and match the hidden layer scale, while the feature linear transformation layer in the back-end of the network is to reduce the feature dimension and match the classification category. The specific calculation of the feature linear transformation layer is as follows (2): x′=xW......(2), Among them, x represents the input feature of a sample, x′ represents the feature of the sample after feature linear transformation, W represents the matrix for feature linear transformation, and W is also the only parameter to be learned in the feature linear transformation layer; 1.3, Frequency domain convolution based on heat diffusion The heat diffusion equation describes the diffusion of heat in a certain space. When this equation is extended to the non-Euclidean domain, the specific equation is as follows (3): Where Δ represents the Laplace-Beltrami operator; u(t) represents the function of the heat at each location in the model space at time t. In the case of discreteness, the LB operator Δ is simplified to the form of an LB matrix, and u(t) is simplified from a function to a finite-dimensional vector. The heat diffusion equation (3) is regarded as a system of first-order linear homogeneous differential equations with constant coefficients. When the boundary condition is t = 0 and u(0) is known, the solution of u(t) is obtained as follows (4): u(t)=e -tΔ u(0)......(4), Among them, e -tΔ The matrix index of e is as follows: After performing eigenvalue decomposition on the matrix A, equation (5) is further simplified. The LB matrix Δ is a real symmetric matrix. Perform eigenvalue decomposition on it. Let the matrix composed of the eigenvectors of the Laplace matrix be U, and the diagonal matrix composed of its corresponding eigenvalues ​​be Λ. Then the Laplace matrix is ​​transformed into Δ = UΛU T , we get the following formula (6): u(t)=Ue -tΛ U T u(0)......(6), Where -tΛ is a diagonal matrix, and the matrix exponentiation of the diagonal matrix e is converted to the element exponentiation directly to simplify the process of solving u(t); Convolution is defined by the heat diffusion equation shown above. Convolution can realize feature fusion in the spatial domain and feature fusion between channels. The above formula (6) is used to realize feature fusion in the spatial domain. The initial feature is set to u(0). The diffusion time t is the parameter to be learned. The size of the fusion area is implicitly determined by t. When t is small, the features in a small area are fused; when t is large, the features in a larger range are fused; when t tends to infinity, the features of all areas in the space are fused. This is similar to the definition of frequency domain convolution. -tΛ It is regarded as a convolution kernel in the frequency domain, and the whole process is considered as a convolution in the frequency domain; 1.

4. Establish a learnable gradient feature module The heat diffusion equation is isotropic, which limits the expressive power of the convolution operator. The gradient of the calculated features is introduced to increase its expressive power and realize anisotropic feature fusion, as follows: In order to calculate the gradient of the function defined on the mesh vertex domain, the function on the mesh vertex domain is regarded as a simplified representation of the function defined on the mesh manifold. For the function defined on the manifold, the gradient at any point is a two-dimensional vector located in its tangent plane. The two-dimensional vector points to the direction where the function value changes the most. For convenience, complex numbers are used to represent the gradient at any point. For the function on the mesh manifold, the gradient is the same everywhere on a face, and the gradients between different faces are different. The gradients on the edges and vertices are approximated by the following method: First, the neighborhood vertices are projected onto the current point tangent plane; Then the function is Taylor expanded at the current point to obtain a set of gradient-related equations; Finally, the least squares method is used to solve the equations. The gradient on the vertex is quickly obtained using the pre-calculated gradient matrix G. After the gradient function of each feature function is calculated, the learnable gradient feature is calculated according to the following formula (7): h = tanh(<z,Az> )......(7), Among them, z i ∈C HS represents the gradient corresponding to all hidden layer features at grid vertex i; A represents a parameter matrix to be learned. The parameter matrix can be a real matrix or a complex matrix. Complex matrices are used in subsequent experiments; <·,·〉 represents the calculation of the dot product of the complex numbers at the corresponding positions of two complex vectors. The dot product here refers to the dot product of the two-dimensional vectors represented by the complex numbers; tanh is an activation function that can ensure that the final feature range is [-1,1]; the final calculation result h i ∈R HS It is the learnable gradient feature at vertex i; The learnable gradient feature is to first perform a learnable linear combination on the gradient of the input feature at each grid vertex, then calculate the dot product between the result of the gradient linear combination and the original gradient, and finally use the activation function to process the result of the dot product and use it as the final learnable gradient feature; Since the solution of the heat diffusion equation is calculated in the frequency domain, in order to simplify and unify the calculation process, it is also necessary to calculate the gradient function of the characteristic function in the frequency domain. To this end, it is only necessary to do some processing on the gradient matrix G, let Where U represents the matrix formed by concatenating the eigenvectors of the LB matrix. Then we only need to By multiplying it with the function transformed to the frequency domain, we can get the gradient function of the function in the spatial domain; 1.5, Multilayer Perceptron Module In order to achieve the fusion of feature information between different channels, a multi-layer perceptron (MLP) is added at the end of each diffusion block. The input of MLP consists of three parts, namely the original input feature x, the gradient feature obtained by the learnable gradient feature module, and the diffusion result returned from the diffusion vector calculated in the frequency domain to the spatial domain.

3. According to the method for segmenting intracranial aneurysms in multimodal medical images based on geometric deep learning as described in claim 1, step 2 is specifically: Step 2.1, processing the reconstructed mesh model to make it a manifold; Step 2.2, calculate the point and surface features of the three-dimensional mesh model. The Euclidean space coordinates, wave kernel signature, maximum principal curvature, minimum principal curvature, average geodesic distance, shape diameter function and other features are calculated on the vertices of the mesh model for subsequent learning, where: Euclidean space coordinates are the basic features of the grid model, but in order to train together with other features, they need to be normalized so that each dimension of their three-dimensional coordinates is in the range of [-1,1]; The wave nuclear signature WKS comes from the theory of spectral analysis and is an advanced shape feature. WKS is derived from the Schrödinger equation and describes the average probability distribution of quantum particles of different energy levels at each vertex on the grid. WKS will not change due to rigid body transformations such as translation and rotation of the grid and is a typical intrinsic feature. The maximum principal curvature refers to the maximum curvature value among all possible curves passing through a point, which describes the direction in which the surface is most sharply curved at that point. The minimum principal curvature refers to the minimum curvature value among all possible curves passing through that point, which describes the direction in which the surface is most gently curved at that point. The average geodesic distance AGD describes the average geodesic distance from a vertex on the grid to all other points on the grid. When AGD is larger, it means that the point is more likely to be at the end of the grid model; when AGD is smaller, the point is more likely to be at the center of the grid model. The shape diameter function SDF defines a scalar value at each point of the three-dimensional grid. This value is obtained by calculating the distance to the opposite side of the shape along the internal direction of the shape at that point and averaging it in the neighborhood of that point. The value of the shape diameter function geometrically represents the thickness of the local area of ​​the shape, thereby describing the overall structure and local characteristics of the blood vessel or aneurysm. Step 2.3, calculate the Laplace-Beltrami matrix on the triangular mesh model. The Laplace operator is a differential operator for functions whose domain is Euclidean space. After being acted on, the function will become the divergence of its own gradient. The Laplace-Beltrami operator is also a differential operator for functions whose domain is a manifold. Similar to the Laplace operator, the function after being acted on will also become the divergence of its own gradient. When facing a mesh manifold, the continuous function on it is simplified by the multidimensional vector of the function value at each vertex of the mesh. The function values ​​at other positions can be obtained by interpolation of the function values ​​of the surrounding vertices. At this time, the LB operator for the function defined on this domain also degenerates into the Laplace-Beltrami matrix. The Cotan algorithm is used to calculate the LB matrix. The result obtained by this algorithm is called the Cotan-Laplace matrix. The calculation method of the Cotan-Laplace matrix when acting on any function is as follows (8): Where Δ represents the Cotan-Laplace matrix, u represents a function defined on the grid manifold, (Δu) i It represents the value of u at vertex i after Δ acts on it, that is, the divergence of the gradient of u at vertex i; represents the set of all neighboring vertices of vertex i; α j and β j represents the diagonal angle between vertex i and adjacent vertex j, u i and u j represents the function value of function u at vertex i and vertex j. The above formula describes the calculation method of the Cotan-Laplacian matrix when it acts on a function defined on a mesh manifold. From it, the specific form of the Cotan-Laplacian matrix can be derived as follows: Among them, Δ i，j represents the value at the i-th row and j-th column of the Cotan-Laplacian matrix; the other symbols are consistent with those in equation (8). Thus, the Cotan-Laplacian matrix is ​​a symmetric matrix. According to the properties of symmetric matrices, the Cotan-Laplacian matrix is ​​subjected to eigenvalue decomposition to obtain its eigenvalues ​​and the corresponding eigenvectors. After obtaining the eigenvalues ​​and eigenvectors of the LB matrix, the frequency domain can be defined on the grid manifold in the same way as the frequency domain in the Euclidean space. For a function whose domain is the Euclidean space, the airspace refers to the Euclidean space where the domain is located, and the frequency domain refers to a functional space formed by the eigenfunction of the Laplace operator as the basis. When it is extended to the function on the grid manifold, the airspace refers to the grid manifold, and the frequency domain refers to the space formed by the eigenvectors of the LB matrix.

4. According to the method for segmenting intracranial aneurysms in multimodal medical images based on geometric deep learning as claimed in claim 1, the multiple loss functions are used in step 3 to reduce the impact of unbalanced data sets, specifically: The multiple loss function consists of three parts, namely, BCE loss function, Focal loss function and Dice loss function. The multiple loss function effectively combines the advantages of the above three loss functions and improves the accuracy of intracranial aneurysm segmentation, as shown in the following formula (10): The BCE loss function effectively quantifies the difference between the predicted segmentation results and the actual segmentation results, providing a robust metric for the classification tasks of blood vessels and aneurysms. in, and l i Represent the predicted label and true label of the i-th sample respectively, and N represents the batch size during training. In addition, the Focal loss function enhances the segmentation performance of intracranial aneurysms by adjusting the weights of simple samples and difficult samples, as well as the weights of the vascular part and the aneurysm part, as shown in the following formula (12): in, l i and N are defined the same as in the BCE loss function. The weighting factor α balances the weights of vessels / aneurysms, while the focusing factor γ controls the balance of easy / difficult examples, where α and γ are specified as 0.8 and 1.5, respectively; Finally, the Dice loss function alleviates the impact of the imbalanced ratio of blood vessels and aneurysms, while enhancing the sensitivity at their segmentation boundaries. The Dice loss function is shown in the following formula (13): in, and l i Same definition as formula (12).

Images