Personalized multi-organ grid reconstruction method based on template deformation

Abstract

The invention provides a personalized multi-organ grid reconstruction method based on template deformation, belongs to the technical field of digital medical treatment, medical image processing, digital twinning and deep learning, and is characterized in that global features and local features of an image can be directly extracted from three-dimensional medical image data of a patient to gradually deform a template grid into a personalized grid. The method has the advantages that manual intervention is not needed, and the personalized multi-organ mesh model (including a surface mesh and a tetrahedral mesh) of the patient can be directly reconstructed from the three-dimensional medical image. The model provided by the invention is an end-to-end architecture, and is easy to train, and meanwhile, the lost design can enable the model to converge quickly. The method has good application prospects in the aspects of individualized organ visualization, virtual operation planning, physical simulation, digital twinning and the like of future patients.

Description

Technical Field

[0001] The present invention belongs to the fields of digital medicine, medical image processing, digital twins, and deep learning technology, and specifically relates to a personalized multi-organ mesh reconstruction method based on template deformation, which mainly utilizes global features and local features extracted from the patient's three-dimensional medical images to iteratively deform the multi-organ template mesh to accurately realize the patient's personalized multi-organ mesh reconstruction. Background Art

[0002] Reconstructing personalized multi-organ meshes from 3D medical images has broad application prospects, such as personalized organ morphology visualization, virtual surgery simulation planning, personalized simulation, etc. However, most of the current methods for obtaining personalized mesh models are divided into two steps. First, the organ labels (Segmentation) are obtained from the patient's 3D medical images, and then the marching cube algorithm is used to generate the patient's personalized surface mesh model. If a tetrahedral mesh is required, the tetrahedral segmentation algorithm is also required to further convert the surface mesh into a tetrahedral mesh model.

[0003] Among them, the methods of obtaining labels from three-dimensional medical images can be roughly divided into manual annotation and (semi-) automatic segmentation. One method relies on radiologists to manually annotate the multi-organ structures in three-dimensional medical images. However, this method is not only time-consuming and laborious, but the quality of the labels obtained is heavily dependent on the professional level of the annotating physician and the physician's subjective judgment. Another method is to automatically segment three-dimensional medical images through segmentation algorithms such as deep learning. This method can greatly shorten the time to obtain personalized segmentation labels for patients, but this type of segmentation algorithm often introduces topological errors and easily produces unrealistic and discontinuous labels. In addition, the accuracy of image labels obtained by manual annotation or automatic segmentation depends on the resolution of the medical image. Therefore, the subsequent generated personalized mesh model often requires the introduction of physicians for manual correction.

[0004] In recent years, with the development of related technologies such as geometric deep learning, it has become possible to directly obtain a patient's personalized mesh model from the patient's personalized medical images. Specifically, a convolutional neural network (CNN) can be used to extract the personalized features of the patient's target organ from three-dimensional medical images such as computed tomography (CT) and magnetic resonance imaging (MRI), and then a graph convolutional neural network (GCN) is used to predict the displacement (i.e., deformation) of each vertex on the target organ mesh template prepared in advance based on the features of the image, thereby directly reconstructing the patient's personalized 3D mesh. This method can not only introduce geometric prior knowledge of the reconstructed organ, but also generate a multi-organ mesh model with better anatomical consistency. Summary of the invention

[0005] The purpose of the present invention is to provide a personalized multi-organ mesh reconstruction method based on template deformation, which can directly reconstruct a personalized multi-organ mesh (including surface mesh and tetrahedral mesh) from three-dimensional medical image data. First, a target organ mesh of a real patient is obtained as a template by manual annotation and the Fourier position coding features of the surface vertices of the mesh template are calculated. The personalized three-dimensional medical image data of different patients (such as CT, MR images, etc.) are input into a convolutional neural network to obtain the personalized global features and local features of the patients. Then, the features at the vertices of the mesh template are obtained by trilinear interpolation according to the coordinates of the vertices on the mesh template surface. The features are spliced ​​with the Fourier position coding features of the vertices on the mesh template surface, and then the obtained template mesh surface vertex features are input into a graph attention convolutional neural network to predict the displacement (i.e., deformation) of the mesh template vertices. Finally, the predicted displacement is added to the mesh template mesh vertex coordinates to obtain the patient's personalized three-dimensional mesh model.

[0006] The technical solution of the present invention is as follows:

[0007] A personalized multi-organ mesh reconstruction method based on template deformation includes the following steps:

[0008] A. Grid template acquisition and related processing

[0009] A1. Generation of multi-organ surface mesh template and multi-organ tetrahedron mesh template: According to the need, manually mark the target organ structure in the template image to obtain the target organ label. Use the marching cube algorithm to obtain the initial multi-organ surface mesh template (You can also add operations such as surface smoothing to make the obtained multi-organ surface mesh template more accurate), then perform tetrahedral mesh division on the initial multi-organ surface mesh template to obtain a multi-organ tetrahedral mesh template Where V sur Represents the vertex coordinates of the multi-organ surface mesh template, E sur represents the edge of the multi-organ surface mesh template, that is, the connection relationship between vertices, n represents the number of vertices in the multi-organ surface mesh template, d represents the number of edges in the multi-organ surface mesh template, V tet represents the vertex coordinates of the multi-organ tetrahedral mesh template, E tet represents the edge of the multi-organ tetrahedron mesh template, m represents the number of vertices of the multi-organ tetrahedron mesh template, and f represents the number of edges in the multi-organ tetrahedron mesh template;

[0010] A2. Calculation of Fourier position coding features of vertices of the multi-organ surface mesh template: Mapping the vertex coordinates of the multi-organ surface mesh template obtained in step A1 to a high-dimensional space by means of random Fourier feature mapping;

[0011] A3. Mesh template deformation based on control point deformation

[0012] A31. Control point acquisition: Use the farthest point sampling algorithm to obtain the multi-organ surface mesh template M from step A1. sur Upsample c vertices and find the multi-organ tetrahedral mesh template M obtained in step A1 with these c vertices tet The coordinates of the nearest point in and serve as control points for mesh deformation of multi-organ surface mesh template and multi-organ tetrahedral mesh template;

[0013] A32, multi-organ tetrahedral mesh template control point deformation weights and multi-organ surface mesh template control point deformation weights Calculation: For the control point vertex coordinates obtained from step A31 Deformation weights of control points of multi-organ tetrahedral mesh template is a linear mapping, that is: V tet =W tet S. Among them, The vertex position V is based on the following tet The optimal solution of the equality constraints with respect to the control point positions yields:

[0014]

[0015] In the formula, is a binary selector matrix used to select the rows of control points in X. is the discrete Bilaplacian matrix of the multi-organ tetrahedral mesh template. When the positions of the selected control points are fixed to S, this optimization can find the vertices V of the multi-organ tetrahedral mesh template. tet The position of the Laplace energy square is minimized. Since the Bilaplacian matrix A is semi-positive definite, this optimization is a convex quadratic programming problem and can be solved as a linear system. As described by Wang et al. in Linear Subspace Design for Real-Time Shape Deformation published in ACM Transactions on Graphics: The solution has V tet =W tet S, where

[0016] W tet =Q T -T T (TAT T ) -1 TAQ T

[0017] In the formula, represents the binary selector matrix complementary to Q.

[0018] Get the tetrahedral mesh control point deformation weight W tet Then, the deformation weights W of the control points of the multi-organ surface mesh template are calculated by interpolation. sur .

[0019] B. Image Data Preprocessing

[0020] B1, data enhancement: including random cropping, random scaling, and random B-spline deformation;

[0021] B2. Image resampling: Method 1: resample the input 3D medical image data voxels to the average value in the three directions of x, y, and z, and then perform center cropping or padding to a fixed size from the resampled image; or; Method 2: resample all 3D medical image data to a fixed voxel size, and then manually or automatically extract the target organ region of interest;

[0022] B3. Image normalization: Normalize the voxel values ​​of the resampled three-dimensional medical image data to between [-1,1].

[0023] C. Reconstruction of personalized multi-organ meshes from 3D medical images

[0024] After obtaining the three-dimensional medical image after image data preprocessing, the input image is firstly feature encoded and decoded using a convolutional neural network to obtain the personalized global features and local features of the patient's three-dimensional image, respectively. The personalized global and local features and the vertex Fourier position encoding features calculated in step A2 are fused to obtain the vertex features of the multi-organ surface mesh template; then the vertex features of the multi-organ surface mesh template are input into the graph attention convolutional neural network to predict the displacement of the vertices of the patient's personalized multi-organ surface mesh model compared to the vertices of the multi-organ surface mesh template. The obtained displacement is added to the vertex coordinates of the multi-organ surface mesh template to obtain the personalized multi-organ surface mesh model predicted by the graph attention convolutional neural network. Then, the personalized multi-organ surface mesh and the personalized multi-organ tetrahedral mesh are calculated using the latest coordinates of the control points on the personalized multi-organ surface mesh model.

[0025] C1. Personalized global and local feature extraction of 3D images

[0026] The convolutional neural network is a U-type neural network, including an encoder and a decoder, which are used to extract personalized global features and local features of patient images, respectively. Among them, the encoder structure consists of K encoding blocks, each encoding block is composed of 1 downsampling layer and W residual blocks, and each residual block is composed of N layers of convolutional layers, and the input and output of each residual block are connected using jump connections. Among them, the downsampling of the encoder is only performed through the downsampling layer in the encoding block, and the channel layer of the corresponding multiple is added while downsampling, that is, if the downsampling factor is 1 / j, the feature channel layer becomes j times the previous layer. The decoder structure consists of K-1 decoding blocks and 1 output layer, and each decoding block is composed of 1 upsampling layer and Y residual blocks. Among them, each decoder enlarges the size of the feature map input of the previous layer by j times through the upsampling layer, and then splices the obtained feature map with the feature map of the corresponding size of the encoder in the channel dimension, and uses the convolution layer to reduce the number of feature channels of the spliced ​​feature map to 1 / j of the number of input channels. At the same time, the decoder layer of the last layer needs to restore the feature map size to the input 3D medical image size and change the number of feature channels to (the current number of organs to be reconstructed + 1) * 2. The output layer uses 1*1*1 convolution to change the number of feature channels to the current number of organs to be reconstructed + 1. The personalized global features of the patient's 3D image refer to the feature map output by the last layer of the encoder, and the local features are the feature maps output by the output layer.

[0027] C2. Personalized multi-organ mesh reconstruction

[0028] C21, Surface vertex feature sampling of multi-organ surface mesh model

[0029] After obtaining the personalized global and local features of the patient's 3D image, the vertex coordinates of the multi-organ surface mesh template obtained in step A1 are used. Perform trilinear interpolation on the global and local features to obtain the personalized features of each vertex in the multi-organ surface mesh template, and then perform feature splicing on the channel dimension with the vertex Fourier position coding features calculated in step A2 to obtain the feature vector of each vertex in the multi-organ surface mesh template. This feature vector is also used as the node feature of the subsequent graph attention convolutional neural network, where n represents the number of vertices of the multi-organ surface mesh template, and F represents the dimension of the feature vector of each vertex of the multi-organ surface mesh template.

[0030] C22, Graph Attention Convolutional Neural Network Predicts Vertex Displacement of Multi-Organ Surface Mesh Templates

[0031] For multi-organ surface mesh template M sur =(V sur ,E sur ), can be represented by a graph structure, where the vertices V of the multi-organ surface mesh template sur is a graph structure node, and the edge of the multi-organ surface mesh template is used as the edge of the graph structure. Therefore, the goal is to predict the displacement of each vertex using the graph attention convolutional neural network based on the feature vector H of each vertex of the multi-organ surface mesh template obtained in step C21

[0032] The graph attention convolutional neural network is composed of k graph attention convolutional layers and 2*(k-1) graph attention residual blocks nested and stacked, where k single-layer graph attention convolutional layers gradually map the feature vector of the vertex from F dimensions to 3 dimensions, that is, the displacement Δs of each vertex in the predicted multi-organ surface mesh template. Each graph attention residual block is composed of y graph attention convolutional layers. It is worth noting that the graph attention convolutional layer in the graph attention residual block does not change the feature dimension. The functions of the graph attention convolutional layer are introduced below:

[0033] In the graph attention convolution layer, the weights of vertices and edges between vertices are shared through a learnable weight Matrix implementation, the weight of the edge between vertices is also called the attention coefficient, the attention coefficient a ij The definition is as follows:

[0034]

[0035] Where σ represents the activation function (such as LeakyReLU, etc.), represents the attention weight parameter learned by the network, and Represent the feature vectors of node i, node j and node b respectively, It represents feature concatenation, N i Represents the set of all vertices connected to vertex i.

[0036] After obtaining the attention coefficient, D independent graph attention mechanisms are further used to enhance the feature learning ability of the network. The graph attention convolution layer outputs the feature for:

[0037]

[0038] Among them, N i represents the set of all vertices connected to vertex i. W is the learnable weight in the graph attention convolutional layer.

[0039] Since the output of the last layer of the graph attention convolutional neural network is the displacement Δs of each vertex in the predicted multi-organ surface mesh template, the output feature of the last graph attention convolutional layer is for:

[0040]

[0041] C23, iterative learning deformation displacement

[0042] After step C22, the graph attention convolutional neural network can predict the displacement Δs of the patient's personalized multi-organ surface mesh relative to the multi-organ surface mesh template, and obtain the vertex coordinates of the personalized multi-organ surface mesh and personalized control point coordinates P c G Then, according to the control point deformation weight W calculated by A32 sur and W tet , the vertex coordinates of the personalized multi-organ surface mesh based on control point deformation can be calculated and personalized multi-organ tetrahedral mesh vertex coordinates

[0043] C24. Iterative Learning

[0044] Vertex coordinates of personalized multi-organ surface mesh obtained using C23 Repeat steps C21 and C22 Z times, and combine the edge E of the multi-organ surface mesh template obtained in step A1 sur and the edge E of the multi-organ tetrahedral mesh template tet , and finally obtain the personalized multi-organ surface mesh model based on control point deformation And the final reconstructed personalized multi-organ tetrahedral mesh model

[0045] C3. Model loss function design

[0046] The weight parameters of the model are trained using the gradient descent algorithm, and the patient-personalized 3D surface mesh model in the training data is used as the ground truth to supervise the model. The loss function is divided into geometric loss, geometric consistency loss, and tetrahedral mesh loss.

[0047] C31, geometric loss

[0048] The geometric loss consists of vertex loss, normal vector loss, edge length loss and Laplace loss.

[0049] (1) Vertex loss

[0050] If the vertex number and topological relationship of the true value label and the multi-organ surface mesh template obtained in step A1 are consistent, the vertex P of the multi-organ surface mesh model of the i-th organ can be predicted by calculation i and the real corresponding organ vertex G i The mean square error between them is taken as the vertex loss L point (P i ,G i ):

[0051]

[0052] If the vertex number or topological relationship of the true value label and the multi-organ surface mesh template obtained in step A1 is inconsistent, the chamfer distance is used as the vertex loss L point (P i ,G i ):

[0053]

[0054] Where p and g represent the vertex predicted by the model and the vertex coordinates in the true label, respectively. In addition, for ease of calculation, the center point of the predicted mesh model patch can be used to replace the predicted mesh model vertex.

[0055] (2) Normal vector loss

[0056] Calculate the normal vector of the predicted multi-organ surface mesh model triangle patch and the normal vector n of the nearest point of the true value label g The loss between:

[0057]

[0058] Where p1 and p2 represent the coordinates of the other two vertices that form a triangular patch with vertex p in the personalized multi-organ surface mesh model predicted by the model.

[0059] (3) Edge length loss

[0060] The average edge length of the predicted personalized multi-organ surface mesh model and the corresponding structure of the ground truth label Loss:

[0061]

[0062] Among them, N p Represents the set of all vertices connected to vertex p.

[0063] (4) Laplace loss

[0064] Calculate the vertex p and its adjacent vertex k p The loss between the average positions is:

[0065]

[0066] In summary, the geometric loss is defined as:

[0067]

[0068] Where P represents the predicted personalized multi-organ surface mesh model, G represents the true value label, N represents the number of organs in the personalized multi-organ surface mesh model, λ1~λ4 represent the weights of vertex loss, normal vector loss, edge length loss, and Laplace loss, respectively, and w i It represents the weight of the loss of different organs.

[0069] C32, geometric consistency loss

[0070] Will and The mean square error between is used as the loss function of geometric consistency loss:

[0071]

[0072] C33, Tetrahedral Mesh Loss

[0073] The tetrahedral mesh loss is defined as follows:

[0074]

[0075] In the formula, Represents the grayscale value of the predicted personalized multi-organ tetrahedral mesh vertex at the corresponding position of the input personalized image, Represents the grayscale value of the personalized multi-organ tetrahedral mesh vertex at the corresponding position in the template image, represents the gray value of the tetrahedron mesh template vertex at the corresponding position of the template image in step A1. Similarly, and They represent the label values ​​of the predicted personalized multi-organ tetrahedral mesh vertices at the corresponding positions of the input personalized image and the label values ​​of the multi-organ tetrahedral mesh template vertices at the corresponding positions of the template labels. α is a hyperparameter of the tetrahedral mesh loss.

[0076] C34, the total loss of the model is:

[0077]

[0078] Where λ5 and λ6 are the weights of the geometric consistency loss and tetrahedral mesh loss respectively, and the weights can be adjusted according to different objectives.

[0079] The effect and benefit of the present invention is that global features and local features can be learned directly from three-dimensional medical images without manual intervention and iterative sampling can be performed, gradually transforming the template grid into a personalized multi-organ grid model of the patient. The model proposed by the present invention is an end-to-end architecture, which is easy to train and the loss design can allow the model to converge quickly. The present invention has a good application prospect in the future for personalized organ visualization of patients, virtual surgery planning, physical simulation, digital twins and other directions. DETAILED DESCRIPTION

[0080] The following takes the patient's personalized heart multi-organ mesh reconstruction as an example, and directly reconstructs a personalized whole heart mesh model from the patient's heart CT image. The reconstructed organs include 7 structures, including the left ventricle, left ventricular myocardium, left atrium, right ventricle, right atrium, aorta, and pulmonary artery. For the convenience of expression, the heart multi-organ mesh reconstruction is referred to as whole heart mesh reconstruction. The specific steps are as follows:

[0081] A personalized multi-organ mesh reconstruction method based on template deformation includes the following steps:

[0082] A. Acquisition of full-heart grid template and related processing

[0083] A1. Generation of whole heart surface mesh template and whole heart tetrahedron mesh template: manually annotate the target organ structure in the template image to obtain labels of seven structures, including left ventricle, left ventricular myocardium, left atrium, right ventricle, right atrium, aorta, and pulmonary artery. Use the marching cube algorithm and use surface smoothing and meshing to obtain the initial whole heart surface mesh template. Then, the obtained full-heart surface mesh template is tetrahedron meshed to obtain the full-heart tetrahedron mesh template Where V sur represents the vertex coordinates of the fully central surface mesh template, E sur represents the edge of the full-heart surface mesh template, that is, the connection relationship between vertices, n represents the number of vertices in the full-heart surface mesh template, d represents the number of edges in the full-heart surface mesh template, V tetrepresents the vertex coordinates of the fully centered tetrahedral mesh template, E tet represents the edge of the fully centered tetrahedron mesh template, m represents the number of vertices of the fully centered tetrahedron mesh template, and f represents the number of edges in the fully centered tetrahedron mesh template;

[0084] A2. Calculation of Fourier position coding features of vertices of the full-heart surface mesh template: Mapping the vertex coordinates of the full-heart surface mesh template obtained in step A1 to a high-dimensional space by using a random Fourier feature mapping method, with the dimension set to 384 dimensions;

[0085] A3. Mesh template deformation based on control point deformation

[0086] A31. Control point acquisition: Use the farthest point sampling algorithm to obtain the full-heart surface mesh template M from step A1. sur Upsample c = 600 vertices and find the multi-organ tetrahedron mesh template M obtained in step A1 with these c vertices tet The coordinates of the nearest point in And it serves as the control point for mesh deformation of the fully-centered surface mesh template and the fully-centered tetrahedral mesh template;

[0087] A32, fully centered tetrahedron mesh template control point deformation weights and the full-heart surface mesh template control point deformation weights Calculation: For the control point vertex coordinates obtained from step A31 Fully centered tetrahedron mesh template control point deformation weights is a linear mapping, that is: V tet =W tet S. Among them, The vertex position V is based on the following tet The optimal solution of the equality constraints with respect to the control point positions yields:

[0088]

[0089] In the formula, is a binary selector matrix used to select the rows of control points in X. is the discrete Bilaplacian matrix of the fully centered tetrahedral mesh template. When the positions of the selected control points are fixed to S, this optimization can find the vertices V of the fully centered tetrahedral mesh template tetThe position of the Laplace energy square is minimized. Since the Bilaplacian matrix A is semi-positive definite, this optimization is a convex quadratic programming problem and can be solved as a linear system. As described by Wang et al. in Linear Subspace Design for Real-Time ShapeDeformation published in ACM Transactions on Graphics: The solution has V tet =W tet S, where

[0090] W tet =Q T -T T (TAT T ) -1 TAQ T

[0091] In the formula, represents the binary selector matrix complementary to Q.

[0092] Get the deformation weight W of the fully centered tetrahedral mesh control points tet Then, the deformation weight W of the control point of the full-heart surface mesh template is calculated by interpolation. sur .

[0093] B. Cardiac CT Image Data Preprocessing

[0094] B1, data enhancement: including random cropping, random scaling, and random B-spline deformation;

[0095] B2. Image resampling: resample the input cardiac CT image data voxels to the average value in the three directions of x, y, and z, and then perform center cropping or padding on the resampled image to a size of [128, 128, 128];

[0096] B3. Image normalization: Normalize the voxel values ​​of the resampled three-dimensional medical image data to between [-1,1].

[0097] C. Reconstruction of personalized whole-heart grid from CT images

[0098] After obtaining the three-dimensional medical image after image data preprocessing, the input image is firstly feature encoded and decoded using a convolutional neural network to obtain the personalized global features and local features of the patient's three-dimensional image, respectively. The personalized global and local features and the vertex Fourier position encoding features calculated in step A2 are fused to obtain the vertex features of the multi-organ surface mesh template; then the vertex features of the multi-organ surface mesh template are input into the graph attention convolutional neural network to predict the displacement of the vertices of the patient's personalized multi-organ surface mesh model compared to the vertices of the multi-organ surface mesh template. The obtained displacement is added to the vertex coordinates of the multi-organ surface mesh template to obtain the personalized multi-organ surface mesh model predicted by the graph attention convolutional neural network. Then, the personalized multi-organ surface mesh and the personalized multi-organ tetrahedral mesh are calculated using the latest coordinates of the control points on the personalized multi-organ surface mesh model.

[0099] C1. Personalized global and local feature extraction of 3D images

[0100] The convolutional neural network is a U-type neural network, including an encoder and a decoder, which are used to extract personalized global features and local features of patient images, respectively. The encoder structure consists of K = 5 encoding blocks, each encoding block is composed of a downsampling layer composed of 3*3*3 convolutions and W = 2 residual blocks, and each residual block is composed of N = 2 layers of 3*3*3 convolution layers, and the input and output of each residual block are connected using a jump connection. The downsampling layer of the encoder also uses 3*3*3 convolutions, each time reducing the input features to the original 1 / j = 1 / 2, and the number of feature channels becomes j = 2 times that of the previous layer. The number of feature map channels output by each encoder is 16, 32, 64, 128 and 256, respectively. The decoder structure consists of K-1 = 4 decoding blocks and 1 output layer, and each decoding block consists of an upsampling layer and a residual block. Among them, each decoder enlarges the size of the feature map input by the previous layer by j = 2 times through the upsampling layer, and then concatenates the obtained feature map with the feature map of the corresponding size of the encoder in the channel dimension, and uses the concatenated feature map to reduce the number of feature channels to 1 / j = 1 / 2 of the number of input channels using the convolution layer. The number of feature channels output by each decoding block is 128, 64, 32 and 16 respectively. At the same time, the decoder layer of the last layer needs to restore the size of the feature map to the input 3D medical image size [128,128,128], and change the number of feature channels to the number of current organs to be reconstructed + 1 = 16. The output layer uses 1*1*1 convolution to change the number of feature channels to the number of target organs to be reconstructed + 1 = 8. The personalized global features of the patient's cardiac CT image refer to the feature map output by the last layer of the encoder (the number of feature map channels is 256), and the local features are the feature maps output by the output layer (the number of channels is 8).

[0101] C2. Personalized whole-heart mesh reconstruction

[0102] C21, full-heart surface mesh model surface vertex feature sampling

[0103] After obtaining the personalized global and local features of the patient's cardiac CT image, the vertex coordinates of the whole heart surface mesh template obtained in step A1 are used. Perform trilinear interpolation on the global and local features to obtain the personalized features of each vertex in the full-heart surface mesh template, and then perform feature splicing on the channel dimension with the vertex Fourier position coding features calculated in step A2 to obtain the feature vector of each vertex in the full-heart surface mesh template. This feature vector is also used as the node feature of the subsequent graph attention convolutional neural network, where n represents the number of vertices of the multi-organ surface mesh template, and F represents the dimension of the feature vector of each vertex of the multi-organ surface mesh template.

[0104] C22, Graph Attention Convolutional Neural Network Predicts Vertex Displacement of Full-Heart Surface Mesh Template

[0105] For the full-heart surface mesh template M sur =(V sur ,E sur ), can be represented by a graph structure, where the vertices V of the multi-organ surface mesh template sur is a graph structure node, and the edge of the full-heart surface mesh template is used as the edge of the graph structure. Therefore, the goal is to predict the displacement of each vertex using the graph attention convolutional neural network based on the feature vector H of each vertex of the full-heart surface mesh template obtained in step C21

[0106] The graph attention convolutional neural network is composed of k=3 graph attention convolutional layers and 2*(k-1)=4 graph attention residual blocks nested and stacked, where k=3 single-layer graph attention convolutional layers gradually map the feature vector of the vertex from F dimensions to 3 dimensions, that is, the displacement Δs of each vertex in the predicted multi-organ surface mesh template. Each graph attention residual block is composed of z=2 graph attention convolutional layers. It is worth noting that the graph attention convolutional layer in the graph attention residual block does not change the feature dimension. The functions of the graph attention convolutional layer are introduced below:

[0107] In the graph attention convolution layer, the weights of vertices and edges between vertices are shared through a learnable weight Matrix implementation, the weight of the edge between vertices is also called the attention coefficient, the attention coefficient a ij The definition is as follows:

[0108]

[0109] Where σ represents the activation function (such as LeakyReLU, etc.), represents the attention weight parameter learned by the network, and Represent the feature vectors of node i, node j and node b respectively, It represents feature concatenation, N i Represents the set of all vertices connected to vertex i.

[0110] After obtaining the attention coefficient, D = 3 independent graph attention mechanisms are further used to enhance the feature learning ability of the network. The graph attention convolution layer outputs the feature for:

[0111]

[0112] Among them, N i represents the set of all vertices connected to vertex i. W is the learnable weight in the graph attention convolutional layer.

[0113] Since the output of the last layer of the graph attention convolutional neural network is the displacement Δs of each vertex in the predicted multi-organ surface mesh template, the output feature of the last graph attention convolutional layer is for:

[0114]

[0115] C23, iterative learning deformation displacement

[0116] After step C22, the graph attention convolutional neural network can predict the displacement Δs of the patient's personalized whole heart surface mesh relative to the whole heart surface mesh template, and obtain the vertex coordinates of the personalized whole heart surface mesh and personalized control point coordinates P c G Then, according to the control point deformation weight W calculated by A32 sur and W tet , the vertex coordinates of the personalized full-heart surface mesh based on the control point deformation can be calculated and personalized fully centered tetrahedron mesh vertex coordinates

[0117] C24. Iterative Learning

[0118] Vertex coordinates of the personalized full-heart surface mesh obtained using C23 Repeat steps C21 and C22 Z = 3 times, and combine the edge E of the full-heart surface mesh template obtained in step A1 sur and the edge E of the fully centered tetrahedral mesh template tet , and finally reconstruct the personalized full-heart surface mesh model based on control point deformation And the final reconstructed personalized fully centered tetrahedron mesh model

[0119] C3. Model loss function design

[0120] The weight parameters of the model are trained using the gradient descent algorithm, and the patient-personalized 3D surface mesh model in the training data is used as the ground truth to supervise the model. The loss function is divided into geometric loss, geometric consistency loss, and tetrahedral mesh loss.

[0121] C31, geometric loss

[0122] The geometric loss consists of vertex loss, normal vector loss, edge length loss and Laplace loss.

[0123] (1) Vertex loss

[0124] Since the vertex number or topological relationship of the true value label in the public dataset and the full-heart surface mesh template obtained in step A1 are inconsistent, this experiment uses the chamfer distance as the vertex loss L point (P i ,G i ):

[0125]

[0126] Where p and g represent the vertex predicted by the model and the vertex coordinates in the true label, respectively. In addition, for ease of calculation, the center point of the predicted mesh model patch can be used to replace the predicted mesh model vertex.

[0127] (2) Normal vector loss

[0128] Calculate the normal vector of the predicted triangular patch of the full-heart surface mesh model and the normal vector n of the nearest point of the true value label g The loss between:

[0129]

[0130] Where p1 and p2 represent the coordinates of the other two vertices that form a triangular patch with vertex p in the personalized full-heart surface mesh model predicted by the model.

[0131] (3) Edge length loss

[0132] The average edge length of the predicted personalized full heart surface mesh model and the corresponding structure of the ground truth label Loss:

[0133]

[0134] Among them, N pRepresents the set of all vertices connected to vertex p.

[0135] (4) Laplace loss

[0136] Calculate the vertex p and its adjacent vertex k p The loss between the average positions is:

[0137]

[0138] In summary, the geometric loss is defined as:

[0139]

[0140] Where P represents the predicted personalized whole heart surface mesh model, G represents the true value label, N represents the number of organs in the personalized whole heart surface mesh model, λ1~λ4 represent the weights of vertex loss, normal vector loss, edge length loss, and Laplace loss, which are 0.3, 0.46, 0.16, and 0.05, respectively, and w i It represents the loss weight of different organs, and all size settings are set to 1.

[0141] C32, geometric consistency loss

[0142] Will and The mean square error between is used as the loss function of geometric consistency loss:

[0143]

[0144] C33, Tetrahedral Mesh Loss

[0145] The tetrahedral mesh loss is defined as follows:

[0146]

[0147] In the formula, Represents the grayscale value of the predicted personalized fully centered tetrahedron mesh vertex at the corresponding position of the input personalized image, Represents the gray value of the personalized fully centered tetrahedron mesh vertex at the corresponding position in the template image, represents the gray value of the tetrahedron mesh template vertex at the corresponding position of the template image in step A1. Similarly, and They represent the label values ​​of the predicted personalized all-centered tetrahedron mesh vertices at the corresponding positions of the input personalized image and the label values ​​of the all-centered tetrahedron mesh template vertices at the corresponding positions of the template labels. α is a hyperparameter of the tetrahedron mesh loss.

[0148] C34, the total loss of the model is:

[0149]

[0150] Where λ5 and λ6 are the weights of the geometric consistency loss and tetrahedral mesh loss, and the weights are set to 10 and 50 respectively.

Claims

1. A personalized multi-organ mesh reconstruction method based on template deformation, characterized in that: The steps include: A. Grid template acquisition and related processing A1. Generation of multi-organ surface mesh template and multi-organ tetrahedron mesh template: According to the needs, the target organ structure is annotated in the template image by manual annotation to obtain the target organ label; Using the marching cubes algorithm to obtain the initial multi-organ surface mesh template Then, the initial multi-organ surface mesh template is tetrahedron-meshed to obtain a multi-organ tetrahedron mesh template. Where V sur Represents the vertex coordinates of the multi-organ surface mesh template, E sur represents the edge of the multi-organ surface mesh template, that is, the connection relationship between vertices, n represents the number of vertices in the multi-organ surface mesh template, d represents the number of edges in the multi-organ surface mesh template, V tet represents the vertex coordinates of the multi-organ tetrahedral mesh template, E tet represents the edge of the multi-organ tetrahedron mesh template, m represents the number of vertices of the multi-organ tetrahedron mesh template, and f represents the number of edges in the multi-organ tetrahedron mesh template; A2. Calculation of Fourier position coding features of vertices of the multi-organ surface mesh template: Mapping the vertex coordinates of the multi-organ surface mesh template obtained in step A1 to a high-dimensional space by means of random Fourier feature mapping; A3. Mesh template deformation based on control point deformation A31. Control point acquisition: Use the farthest point sampling algorithm to obtain the multi-organ surface mesh template M from step A1. sur Upsample c vertices and find the multi-organ tetrahedral mesh template M obtained in step A1 with these c vertices tet The coordinates of the nearest point in and serve as control points for mesh deformation of multi-organ surface mesh template and multi-organ tetrahedral mesh template; A32, multi-organ tetrahedral mesh template control point deformation weights and multi-organ surface mesh template control point deformation weights Calculation: For the control point vertex coordinates obtained from step A31 Deformation weights of control points of multi-organ tetrahedral mesh template is a linear mapping, that is: V tet =W tet S; among them, The vertex position V is based on the following tet The optimal solution of the equality constraints with respect to the control point positions yields: In the formula, is a binary selector matrix used to select the rows of control points in X. is the discrete Bilaplacian matrix of the multi-organ tetrahedral mesh template; when the positions of the selected control points are fixed to S, this optimization can find the vertices V of the multi-organ tetrahedral mesh template tet The position minimizes the square of the Laplace energy. Since the Bilaplacian matrix A is semi-positive definite, this optimization is a convex quadratic programming problem and can be solved as a linear system. The solution has V tet =W tet S, where W tet =Q T -T T (TAT T ) -1 TAQ T In the formula, represents a binary selector matrix complementary to Q; Get the tetrahedral mesh control point deformation weight W tet Then, the deformation weights W of the control points of the multi-organ surface mesh template are calculated by interpolation. sur ; B. Image Data Preprocessing B1, data enhancement: including random cropping, random scaling, and random B-spline deformation; B2. Image resampling: Method 1: resample the input 3D medical image data voxels to the average value in the three directions of x, y, and z, and then perform center cropping or padding to a fixed size from the resampled image; or; Method 2: resample all 3D medical image data to a fixed voxel size, and then manually or automatically extract the target organ region of interest; B3, Image normalization: normalize the voxel values ​​of the resampled 3D medical image data to between [-1,1]; C. Reconstruction of personalized multi-organ meshes from 3D medical images After obtaining the three-dimensional medical image after image data preprocessing, first use the convolutional neural network to perform feature encoding and feature decoding on the input image to obtain personalized global features and local features of the patient's three-dimensional image respectively; fuse the personalized global and local features and the vertex Fourier position coding features calculated in step A2 to obtain the vertex features of the multi-organ surface mesh template; then input the vertex features of the multi-organ surface mesh template into the graph attention convolutional neural network to predict the displacement of the vertices of the patient's personalized multi-organ surface mesh model compared to the vertices of the multi-organ surface mesh template; add the obtained displacement to the vertex coordinates of the multi-organ surface mesh template to obtain the personalized multi-organ surface mesh model predicted by the graph attention convolutional neural network; then use the latest coordinates of the control points on the personalized multi-organ surface mesh model to calculate the personalized multi-organ surface mesh and the personalized multi-organ tetrahedral mesh; C1. Personalized global and local feature extraction of 3D images The convolutional neural network is a U-type neural network, including an encoder and a decoder, which are respectively used to extract personalized global features and local features of patient images; wherein the encoder structure is composed of K encoding blocks, each encoding block is composed of 1 downsampling layer and W residual blocks, and each residual block is composed of N convolutional layers, and the input and output of each residual block are connected using jump connections; wherein the encoder downsampling is only performed through the downsampling layer in the encoding block, and the corresponding multiple channel layer is added while downsampling, that is, if the downsampling factor is 1 / j, the feature channel layer becomes j times the previous layer; the decoder structure is composed of K-1 decoding blocks and 1 output layer, each decoding block is composed of 1 upsampling layer and Y residual The decoder is composed of blocks; each decoder enlarges the size of the feature map input by the previous layer by j times through the upsampling layer, and then splices the obtained feature map with the feature map of the corresponding size of the encoder in the channel dimension, and uses the convolution layer to reduce the number of feature channels of the spliced ​​feature map to 1 / j of the number of input channels; at the same time, the decoder layer of the last layer needs to restore the size of the feature map to the size of the input three-dimensional medical image, and change the number of feature channels to (the current number of organs to be reconstructed + 1)*2; the output layer uses 1*1*1 convolution to change the number of feature channels to the current number of organs to be reconstructed + 1. The personalized global features of the patient's three-dimensional image refer to the feature maps output by the last layer of the encoder, and the local features are the feature maps output by the output layer; C2. Personalized multi-organ mesh reconstruction C21, Surface vertex feature sampling of multi-organ surface mesh model After obtaining the personalized global and local features of the patient's 3D image, the vertex coordinates of the multi-organ surface mesh template obtained in step A1 are used. Perform trilinear interpolation on the global and local features to obtain the personalized features of each vertex in the multi-organ surface mesh template, and then perform feature splicing on the channel dimension with the vertex Fourier position coding features calculated in step A2 to obtain the feature vector of each vertex in the multi-organ surface mesh template. This feature vector is also used as the node feature of the subsequent graph attention convolutional neural network, where n represents the number of vertices of the multi-organ surface mesh template, and F represents the dimension of the feature vector of each vertex of the multi-organ surface mesh template; C22, Graph Attention Convolutional Neural Network Predicts Vertex Displacement of Multi-Organ Surface Mesh Templates For multi-organ surface mesh template M sur =(V sur ,E sur ), can be represented by a graph structure, where the vertices V of the multi-organ surface mesh template sur is a graph structure node, and the edge of the multi-organ surface mesh template is used as the edge of the graph structure; therefore, the goal is to predict the displacement of each vertex using the graph attention convolutional neural network based on the feature vector H of each vertex of the multi-organ surface mesh template obtained in step C21 The graph attention convolutional neural network is composed of k graph attention convolutional layers and 2*(k-1) graph attention residual blocks nested and stacked, wherein the k single-layer graph attention convolutional layers gradually map the feature vector of the vertex from F dimension to 3 dimension, that is, the displacement Δs of each vertex in the predicted multi-organ surface mesh template; wherein each graph attention residual block is composed of y graph attention convolutional layers, and it is worth noting that the graph attention convolutional layer in the graph attention residual block does not change the feature dimension; the function of the graph attention convolutional layer is introduced below: In the graph attention convolution layer, the weights of vertices and edges between vertices are shared through a learnable weight Matrix implementation, the weight of the edge between vertices is also called the attention coefficient, the attention coefficient a ij The definition is as follows: In the formula, σ represents the activation function, represents the attention weight parameter learned by the network, and Represent the feature vectors of node i, node j and node b respectively, It represents feature concatenation, N i represents the set of all vertices connected to vertex i; After obtaining the attention coefficient, D independent graph attention mechanisms are further used to enhance the feature learning ability of the network. The graph attention convolution layer outputs the feature for: Among them, N i represents the set of all vertices connected to vertex i; W is the learnable weight in the graph attention convolution layer; Since the output of the last layer of the graph attention convolutional neural network is the displacement Δs of each vertex in the predicted multi-organ surface mesh template, the output feature of the last layer of the graph attention convolutional layer is for: C23, iterative learning deformation displacement After step C22, the graph attention convolutional neural network can predict the displacement Δs of the patient's personalized multi-organ surface mesh relative to the multi-organ surface mesh template, and obtain the vertex coordinates of the personalized multi-organ surface mesh and personalized control point coordinates Then, according to the control point deformation weight W calculated by A32 sur and W tet , the vertex coordinates of the personalized multi-organ surface mesh based on control point deformation can be calculated and personalized multi-organ tetrahedral mesh vertex coordinates C24. Iterative Learning Vertex coordinates of personalized multi-organ surface mesh obtained using C23 Repeat steps C21 and C22 Z times, and combine the edge E of the multi-organ surface mesh template obtained in step A1 sur and the edge E of the multi-organ tetrahedral mesh template tet , and finally obtain the personalized multi-organ surface mesh model based on control point deformation And the final reconstructed personalized multi-organ tetrahedral mesh model C3. Model loss function design The weight parameters of the model are trained using the gradient descent algorithm, and the patient-personalized 3D surface mesh model in the training data is used as the true value label to supervise the model. The loss function is divided into geometric loss, geometric consistency loss, and tetrahedral mesh loss. C31, geometric loss The geometric loss consists of vertex loss, normal vector loss, edge length loss and Laplace loss; (1) Vertex loss If the vertex number and topological relationship of the true value label and the multi-organ surface mesh template obtained in step A1 are consistent, the vertex P of the multi-organ surface mesh model of the i-th organ can be predicted by calculation i and the real corresponding organ vertex G i The mean square error between them is taken as the vertex loss L point (P i ,G i ): If the vertex number or topological relationship of the true value label and the multi-organ surface mesh template obtained in step A1 is inconsistent, the chamfer distance is used as the vertex loss L point (P i ,G i ): Where p and g represent the vertex predicted by the model and the vertex coordinates in the real label respectively; in addition, for the convenience of calculation, the center point of the predicted mesh model patch can be used to replace the predicted mesh model vertex; (2) Normal vector loss Calculate the normal vector of the predicted multi-organ surface mesh model triangle patch and the normal vector n of the nearest point of the true value label g The loss between: Where p1 and p2 represent the coordinates of the other two vertices that form a triangular patch with vertex p in the personalized multi-organ surface mesh model predicted by the model; (3) Edge length loss The average edge length of the predicted personalized multi-organ surface mesh model and the corresponding structure of the ground truth label Loss: Among them, N p represents the set of all vertices connected to vertex p; (4) Laplace loss Calculate the vertex p and its adjacent vertex k p The loss between the average positions is: In summary, the geometric loss is defined as: Where P represents the predicted personalized multi-organ surface mesh model, G represents the true value label, N represents the number of organs in the personalized multi-organ surface mesh model, λ1~λ4 represent the weights of vertex loss, normal vector loss, edge length loss, and Laplace loss, respectively, and w i It represents the weight of the loss of different organs; C32, geometric consistency loss Will and The mean square error between is used as the loss function of geometric consistency loss: C33, Tetrahedral Mesh Loss The tetrahedral mesh loss is defined as follows: In the formula, Represents the grayscale value of the predicted personalized multi-organ tetrahedral mesh vertex at the corresponding position of the input personalized image, Represents the grayscale value of the personalized multi-organ tetrahedral mesh vertex at the corresponding position in the template image, represents the gray value of the tetrahedron mesh template vertex at the corresponding position of the template image in step A1; similarly, and They represent the label values ​​of the predicted personalized multi-organ tetrahedral mesh vertices at the corresponding positions of the input personalized image and the label values ​​of the multi-organ tetrahedral mesh template vertices at the corresponding positions of the template labels respectively; α is the hyperparameter of the tetrahedral mesh loss; C34, the total loss of the model is: Where λ5 and λ6 are the weights of the geometric consistency loss and tetrahedral mesh loss respectively, and the weights can be adjusted according to different objectives.

Images