Neural fiber bundle clustering segmentation method and system based on dispersion supervoxels

Abstract

The invention discloses a nerve fiber bundle clustering segmentation method and system based on diffusion supervoxels, and relates to the technical field of medical image analysis and artificial intelligence, and the method comprises the steps: receiving three-dimensional brain magnetic resonance data, and carrying out the preprocessing; establishing a personalized diffusion magnetic resonance brain template matched with the crowd by using an iterative registration method; establishing a diffusion weighted undirected graph, and realizing diffusion super voxel segmentation based on the diffusion geodesic distance; realizing fiber bundle clustering based on the maximum group of the nerve fiber connection diagram of the dispersion supervoxel; the dispersion super voxels are registered to an individual space, and personalized segmentation of the fiber bundles is realized by using the maximum clique; constructing a center curve of the fiber bundle cluster and mapping the dispersion characteristics to the curve; and establishing a statistical analysis method of a center curve based on replacement test, and establishing a machine learning classification model.

Description

Technical Field

[0001] The present invention relates to the field of medical image analysis and artificial intelligence technology, and in particular to a method and system for clustering and segmenting nerve fiber bundles based on diffuse supervoxels. Background Art

[0002] Diffusion Magnetic Resonance Imaging (dMRI) uses a gradient field to change the direction of free diffusion of water molecules, thereby measuring the degree of free diffusion of water molecules in tissues. Common sequences include Diffusion Tensor Imaging (DTI), Diffusion Kurtosis Imaging (DKI), Diffusion Spectrum Imaging (DSI) and High Angular Resolution Diffusion Imaging (HARDI). The fiber tract tracing method establishes a diffusion tensor model or orientation distribution function (ODF) of water molecule diffusion from dMRI, and then reconstructs the fiber curve that describes the course and direction of white matter fiber tracts.

[0003] The results of fiber tract tracing can reflect the damage of the disease to the white matter fiber tracts, thereby improving the accuracy of diagnosis. Since there are multiple fiber tracts in the brain connecting different brain regions, it is necessary to cluster the fiber tract tracing results to establish a white matter fiber tract atlas, and then segment the individual white matter fiber tracts based on the fiber tract atlas. Current fiber tract clustering methods are mostly based on the distance between fiber curves, such as Hausdorff distance and average shortest distance, which have high computational complexity and require downsampling of fiber tract tracing results. They are not suitable for the construction of fiber tract atlases based on dMRI imaging data of large sample populations, which limits the development of related fields. Summary of the invention

[0004] To address the deficiencies mentioned in the above background technology, the purpose of the present invention is to provide a neural fiber bundle clustering segmentation method and system based on diffuse supervoxels, which converts the white matter fiber bundle clustering problem into a diffuse supervoxel clustering problem, thereby reducing the computational complexity and improving the computational speed.

[0005] In a first aspect, the purpose of the present invention can be achieved by the following technical solution: a neural fiber bundle clustering segmentation method based on diffuse supervoxels, the method comprising the following steps:

[0006] receiving three-dimensional brain magnetic resonance imaging data, performing preprocessing to obtain processed image data, and registering the processed image data using an iterative registration method to obtain a brain template;

[0007] The brain template is input into a pre-established water molecule diffusion model, and a diffusion weighted undirected graph is outputted. The diffusion geodesic distance is generated according to the diffusion weighted undirected graph, and the diffusion supervoxel is obtained based on the diffusion geodesic distance segmentation;

[0008] Taking diffuse supervoxels as vertices, a weighted undirected graph of supervoxel nerve fiber connections is generated, the maximum cluster of the weighted undirected graph of supervoxel nerve fiber connections is extracted, and the weighted undirected graph of supervoxel nerve fiber connections is updated according to the maximum cluster to achieve fiber bundle clustering;

[0009] The diffusion supervoxel is registered to the individual space, and the individual fiber bundles are obtained according to the maximum cluster segmentation of the individual diffusion supervoxel. The fiber density image of the individual fiber bundle is used for curve fitting to obtain the fiber bundle center curve. The characteristics of the fiber bundle center curve are extracted and statistically analyzed to achieve segmentation and analysis of the fiber bundle.

[0010] In combination with the first aspect, in some implementations of the first aspect, the method further includes: the process of registering the processed image data using the iterative registration method:

[0011] Given the processed image data (J1, J2, ..., J n ), the deformations registered to the brain template φ are θ1, θ2, …, θ n , the brain template θ is established by the group registration method:

[0012] The individual images J1, J2, ..., J n Register to current brain template The deformations are t represents the tth iteration;

[0013] Pick The non-rigid deformation part of the inverse transformation is averaged as the residual deformation R (t) ;

[0014] Use residual deformation R (t) Correction deformation in is the corrected deformation;

[0015] Use the corrected deformation Image J i Transform and take the mean of the transformed image as the brain template for the next iteration n is the number of images;

[0016] Repeat the above until the difference between two consecutive generated brain templates is less than the preset value ε:

[0017]

[0018] The brain template construction process adopts a coarse-to-fine approach. First, the brain template is established using affine transformation, and then nonlinear deformation is used to successively establish low-resolution brain templates and high-resolution brain templates.

[0019] Use the brain template to build the process, using the ICBM-152 standard brain template as the initial template Establish T1WI brain template;

[0020] The final step of the T1WI brain template construction process is to correct the deformation The aligned dMRI images are deformed into the T1WI brain template space, and the mean is taken as the initial brain template of dMRI;

[0021] Based on the initial brain template of dMRI, the dMRI brain template is established by applying the above brain template construction process.

[0022] In combination with the first aspect, in certain implementations of the first aspect, the method further includes: a process of inputting the brain template into a pre-established water molecule diffusion model and outputting a diffusion weighted undirected graph:

[0023] Taking white matter voxels as the vertices of the graph, when two voxels have a 26-adjacent relationship, there is an edge e between the corresponding vertices v1 and v2, and the weight of edge e is defined as:

[0024] W e =d diffusion +γd chamfer

[0025] where d diffusion Represents the distance between two voxel points of the water molecule diffusion model, d chamfer Indicates the space chamfer distance;

[0026] For the diffusion tensor model of water molecule diffusion, the difference is defined as the logarithmic Euclidean distance of the diffusion tensor. Assuming that the diffusion tensors corresponding to v1 and v2 are D1 and D2 respectively, the logarithmic Euclidean distance is:

[0027] d tensor =tr((logD1-logD2) T (logD1-logD2))

[0028] Where tr represents the trace of the matrix;

[0029] For the directional distribution function model of water molecule diffusion, the ODF model f can represent a linear combination of spherical harmonic functions:

[0030]

[0031] where θ and are the polar angle and azimuth of the spherical coordinate system, Y lm is the spherical harmonic function, f lm are the spherical harmonic coefficients, l and m are the degree and order of the spherical harmonic function respectively;

[0032] Assume that the ODF models of v1 and v2 are f and g respectively, then according to the Bregman divergence, the square of the L2 norm of the two is:

[0033]

[0034] Where S represents the unit sphere, Ω represents the spherical angle and θ and are the polar angle and azimuth angle of the spherical coordinate system respectively; since the spherical harmonic function is a set of orthogonal bases, the L2 norm between ODFs is simplified to:

[0035]

[0036] In addition, according to the Bregman divergence, the Kullback-Leibler divergence between ODFs is defined as:

[0037]

[0038] Then its symmetric non-negative Jensen-Shannon divergence is d KL (f,g) and d KL The mean of (g,f) is:

[0039]

[0040] Assume lnf lm and lng lm are the spherical harmonic coefficients of lnf and lng respectively. According to the orthogonality of spherical harmonic functions, the Jensen-Shannon divergence is simplified to:

[0041]

[0042] The coefficient γ is used to balance the water molecule diffusion model distance and the space chamfer distance, and is set to γ1 multiplied by the median of the ratio of the water molecule diffusion model distance to the space chamfer distance of the adjacent voxel points.

[0043] In combination with the first aspect, in some implementations of the first aspect, the method further includes: generating a diffuse geodesic distance according to the diffuse weighted undirected graph, initializing a seed point using a KMeans++ algorithm according to the geodesic distance between vertices in the diffuse weighted undirected graph, and the process is as follows:

[0044] Initialization of seed points based on white matter regions is based on diffusion weighted undirected graph G D Select the seed point;

[0045] Based on the seed point initialization of the white matter skeleton, the various fractional anisotropies or generalized fractional anisotropy of the white matter are refined by graphics to obtain the white matter skeleton, and then the G corresponding to the white matter skeleton is extracted. D The subgraph G′ D Select the seed point;

[0046] Assignment step: Calculate the diffuse weighted undirected graph G D The geodesic distance between the vertex and the seed point, assigning the vertex to the cluster corresponding to the seed point with the minimum geodesic distance;

[0047] Seed point update step: Use the dynamic programming algorithm to select the vertex with the smallest sum of geodesic distances to other vertices in the cluster as the new seed point;

[0048] Iterate the assignment step and the seed point update step until the seed point does not change or cycles in a certain set of vertices; the voxels corresponding to the vertices in the same cluster form a diffuse supervoxel.

[0049] In combination with the first aspect, in some implementations of the first aspect, the method further includes: the fiber bundle clustering is performed based on diffuse supervoxels, and the process is as follows:

[0050] Construct a weighted undirected graph G of supervoxel nerve fiber connections with diffuse supervoxels as vertices S , if and only if there exists a supervoxel S i and S j When the nerve fiber bundle is connected, there is an edge e(S i ,S j );

[0051] Assume that after S i and S j The fiber bundle sets are T i and T j , the volumes are V i and V j , then the edge e(S i ,S j ) is:

[0052]

[0053] For Figure G SThe adjacency matrix of is binarized, and the Otsu algorithm can be used to find the threshold, or the sparsity of the adjacency matrix after binarization |E| / K(K-1) satisfies the given conditions, where |E| is the number of edges after binarization, and K is the number of diffuse supervoxels;

[0054] Use Bron-Kerbosch to extract all the largest groups of the binarized graph, which are all the largest connected components;

[0055] Take the largest clique C with the largest internal connection max ;

[0056] Take the least number of passes through C max The fiber bundles with diffuse supervoxels within ξ% are regarded as clusters, which are deleted from all fiber bundles and the graph G is updated. S The weight of the edge, ξ is the preset ratio;

[0057] Iterate the above steps until the newly obtained C max The internal connection is less than a given threshold, where the threshold is K×G S Adjacency matrix element mean.

[0058] In combination with the first aspect, in certain implementations of the first aspect, the method further includes: performing curve fitting according to the fiber density image of the individual fiber bundle to obtain the fiber bundle center curve:

[0059] Up-sample the dMRI brain template and construct fiber bundle density maps of fiber bundle clusters;

[0060] The central curve of the fiber bundle density map is fitted using a non-uniform rational B-spline curve, and the diffusion characteristics of the fiber bundle are mapped to the central curve to obtain the central curve of the fiber bundle.

[0061] In combination with the first aspect, in some implementations of the first aspect, the method further includes: a calculation process of the fiber bundle center curve:

[0062] For any cluster of nerve fiber bundles, the nerve fiber density map was calculated in the upsampled dMRI brain template space, which was the ratio of the number of fiber curves passing through each voxel to the cluster size;

[0063] Use non-uniform rational B-spline curve to fit the nerve fiber density map:

[0064]

[0065] Where P i is the control point, N i,p (u) is the B-spline function, u takes values ​​between 0 and 1, Curve(0) and Curve(1) are the center points of the diffuse supervoxels at the two ends of the largest cluster corresponding to the fiber bundle;

[0066] Let Q0, Q1, …, Q m are the voxel points through which the fiber bundles pass, and their densities are d0, d1, …, d m , NURBS curve fitting uses the least squares method to minimize the objective function:

[0067]

[0068] The NURBS curve is deformed to the individual space using the registered quasi-deformation, and then the diffusion characteristics such as the weighted average FA and GFA of the fiber bundle are mapped to the NURBS curve;

[0069] For the points on the curve, the fiber density map is divided into a Voronoi diagram based on the Euclidean distance, and the diffuse features in the Voronoi diagram are mapped to the corresponding points:

[0070]

[0071] where v k Represents the voxel Q in the Voronoi diagram k diffusion characteristics.

[0072] In a second aspect, in order to achieve the above-mentioned object, the present invention discloses a neural fiber bundle clustering segmentation system based on diffuse supervoxels, comprising:

[0073] An image acquisition and preprocessing module is used to receive three-dimensional brain magnetic resonance image data, perform preprocessing, obtain processed image data, and register the processed image data using an iterative registration method to obtain a brain template;

[0074] A diffusion supervoxel segmentation module is used to input the brain template into a pre-established water molecule diffusion model, output a diffusion weighted undirected graph, generate a diffusion geodesic distance based on the diffusion weighted undirected graph, and obtain a diffusion supervoxel based on the diffusion geodesic distance segmentation;

[0075] A fiber bundle clustering module is used to generate a weighted undirected graph of supervoxel nerve fiber connections with diffuse supervoxels as vertices, extract the maximum cluster of the weighted undirected graph of supervoxel nerve fiber connections, and update the weighted undirected graph of supervoxel nerve fiber connections according to the maximum cluster to achieve fiber bundle clustering;

[0076] The fiber bundle center curve mapping module is used to align the diffusion supervoxel to the individual space, obtain the individual fiber bundle according to the maximum cluster segmentation of the individual diffusion supervoxel, perform curve fitting based on the fiber density image of the individual fiber bundle to obtain the fiber bundle center curve, extract the characteristics of the fiber bundle center curve and perform statistical analysis to realize the segmentation analysis of the fiber bundle.

[0077] In another aspect of the present invention, in order to achieve the above-mentioned purpose, a terminal device is disclosed, including a memory, a processor, and a computer program stored in the memory and capable of running on the processor, wherein the memory stores a computer program capable of running on the processor, and when the processor loads and executes the computer program, the above-mentioned diffuse supervoxel-based neural fiber bundle clustering segmentation method is adopted.

[0078] In another aspect of the present invention, in order to achieve the above-mentioned purpose, a computer-readable storage medium is disclosed, in which a computer program is stored, characterized in that when the computer program is loaded and executed by a processor, the above-mentioned diffuse supervoxel-based neural fiber bundle clustering segmentation method is adopted.

[0079] Beneficial effects of the present invention:

[0080] The present invention converts the complex fiber bundle clustering problem into an easy-to-solve diffuse supervoxel clustering problem, which has the advantages of fast computing speed and less memory required, and can process large-scale fiber curves, providing a feasible solution for the construction of fiber bundle maps for large sample populations. It segments fiber bundles according to the spatial relationship through diffuse supervoxels, and has the advantages of fast computing speed and high robustness, providing a new solution for statistical analysis and classification model training of specific fiber bundles. BRIEF DESCRIPTION OF THE DRAWINGS

[0081] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, for those skilled in the art, other drawings can be obtained based on these drawings without creative work.

[0082] Figure 1 is a flow chart of the method of the present invention;

[0083] Figure 2 is a flow chart of constructing a dMRI brain template in an embodiment of the present invention;

[0084] Figure 3 is a structural block diagram of diffuse supervoxel segmentation in an embodiment of the present invention;

[0085] Figure 4 is a flow chart of fiber bundle clustering based on diffuse supervoxels in an embodiment of the present invention;

[0086] Figure 5 It is a schematic diagram of the system structure of the present invention. DETAILED DESCRIPTION

[0087] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0088] Embodiment 1:

[0089] The following is an introduction to the relevant terms involved in the embodiments of the present application:

[0090] Nerve fiber bundle: In the white matter, nerve fibers with basically the same starting and ending points, course and function are grouped together to form a fiber bundle.

[0091] like Figure 1 As shown, the neural fiber bundle clustering segmentation method based on diffuse supervoxel includes the following steps:

[0092] S101: receiving three-dimensional brain magnetic resonance imaging data, and performing preprocessing to obtain processed image data, and registering the processed image data using an iterative registration method to obtain a brain template;

[0093] The preprocessing steps include removal of non-brain tissue, gray-white matter segmentation, registration, water molecule diffusion model reconstruction, and fiber bundle tracing;

[0094] The process of registering the processed image data using an iterative registration method:

[0095] like Figure 2 As shown, given a set of preprocessed image data (J1, J2, ..., J n ), the deformations registered to the brain template φ are θ1, θ2, …, θ n , the brain template θ can be established by the group registration method:

[0096] i. The individual images J1, J2, ..., J n Register to current brain template The deformations are

[0097] ii. Take The non-rigid deformation part of the inverse transformation is averaged as the residual deformation R (t) ;

[0098] iii. Use residual deformation R (t) Correction deformation

[0099] iv. Use the corrected deformation Image J iTransform and take the mean of the transformed image as the brain template for the next iteration

[0100] v. Repeat steps i–iv until the difference between two consecutive generated brain templates is less than the preset value:

[0101] The brain template construction process adopts a coarse-to-fine approach. First, the brain template is established using affine transformation, and then nonlinear deformation is used to successively establish low-resolution brain templates and high-resolution brain templates.

[0102] Use the above brain template construction process and take the ICBM-152 standard brain template as the initial template Establish T1WI brain template;

[0103] The final step of the T1WI brain template construction process is to correct the deformation The aligned dMRI images are deformed into the T1WI brain template space, and the mean is taken as the initial brain template of dMRI;

[0104] Based on the initial brain template of dMRI, the dMRI brain template is established by applying the above-mentioned brain template construction process.

[0105] S102: inputting the brain template into a pre-established water molecule diffusion model, outputting a diffusion weighted undirected graph, generating a diffusion geodesic distance according to the diffusion weighted undirected graph, and obtaining diffusion supervoxels based on the diffusion geodesic distance segmentation;

[0106] like Figure 3 As shown in Figure 2, diffusion supervoxel segmentation establishes a diffusion weighted undirected graph G according to the difference between the spatial adjacency relationship of voxel points in the white matter brain region and the water molecule diffusion model. D Then, the seed point (i.e., the central voxel of the cluster) is initialized in the white matter region or the white matter skeleton region using the KMeans++ method based on diffusion distance, and finally the white matter region is over-segmented into diffuse supervoxels using a clustering method based on diffusion geodesic distance.

[0107] Diffuse weighted undirected graph G D Build, specifically including:

[0108] Taking white matter voxels as the vertices of the graph, when two voxels have a 26-adjacent relationship (optionally 6-adjacent or 18-adjacent) there is an edge e between their corresponding vertices v1 and v2. The weight of edge e is defined as:

[0109] W e =d diffusion +γd chamfer

[0110] where d diffusionRepresents the distance between two voxel points of the water molecule diffusion model, d chamfer Indicates the space chamfer distance;

[0111] For the diffusion tensor model of water molecule diffusion, the difference is defined as the logarithmic Euclidean distance of the diffusion tensor. Assuming that the diffusion tensors corresponding to v1 and v2 are D1 and D2 respectively, their logarithmic Euclidean distance is:

[0112] d tensor =tr((logD1-logD2) T (logD1-logD2))

[0113] Where tr represents the trace of the matrix;

[0114] For the Orientation Distribution Function (ODF) model of water molecule diffusion, the ODF model f can represent a linear combination of spherical harmonic functions:

[0115]

[0116] where θ and are the polar angle and azimuth of the spherical coordinate system, Y lm is the spherical harmonic function, f lm is the spherical harmonic coefficient, l is the order of the spherical harmonic function;

[0117] Assume that the ODF models of v1 and v2 are f and g respectively, then according to the Bregman divergence, the square of the L2 norm of the two is:

[0118]

[0119] Where S represents the unit sphere, Ω represents the spherical angle and Since spherical harmonics are a set of orthogonal bases, the L2 norm between ODFs is simplified to:

[0120]

[0121] In addition, according to the Bregman divergence, the Kullback-Leibler divergence between ODFs is defined as:

[0122]

[0123] Then its symmetric non-negative Jensen-Shannon divergence is d KL (f,g) and d KL The mean of (g,f) is:

[0124]

[0125] Assume lnflm and lng lm are the spherical harmonic coefficients of lnf and lng respectively. According to the orthogonality of spherical harmonic functions, the Jensen-Shannon divergence can be simplified to:

[0126]

[0127] To simplify the calculation, the upper limit of the order l of the spherical harmonic function can be selected as 4 or 6;

[0128] The coefficient γ is used to balance the water molecule diffusion model distance and the space chamfer distance. It is set to γ1 multiplied by the median of the ratio of the water molecule diffusion model distance to the space chamfer distance of the adjacent voxel point. γ1 can be a real number (0,1].

[0129] Generate diffuse geodesic distance based on diffuse weighted undirected graph Based on the geodesic distance between vertices in the diffuse weighted undirected graph, use KMeans++ algorithm to initialize seed points. The process is as follows:

[0130] Initialization of seed points based on white matter regions is based on diffusion weighted undirected graph G D Select the seed point;

[0131] Based on the seed point initialization of the white matter skeleton, the fractional anisotropy (FA) or generalized fractional anisotropy (GFA) of the white matter is refined by graphics to obtain the white matter skeleton, and then the G corresponding to the white matter skeleton is extracted. D The subgraph G′ D Select a seed point.

[0132] For a given diffuse weighted undirected graph or its subgraph G, the KMeans++ method based on diffuse geodesic distance is used to select seed points:

[0133] i. Initialize G, assign all vertex weights to floating-point maximum values, and mark them as unvisited;

[0134] ii. Randomly select a vertex of G as the first seed point, assign it a weight of 0, mark it as visited and add it to the priority queue Q;

[0135] iii. Pop the head vertex v of the queue Q and visit the unvisited adjacent vertices u of v in turn;

[0136] iv. If w v +w e (v,u) <w u , then mark u as visited and update its weight to w v +w e(v,u) and join the tail of queue Q, where w v and w u are the weights of vertices v and u, respectively, and w e (v,u)

[0137] is the weight of the edge connecting vertices v and u;

[0138] v. Repeat steps iii and iv until the queue Q is empty, and the vertex weight is the diffuse geodesic distance from the nearest seed point;

[0139] vi. For all non-seed vertices, establish a distribution based on the square of the vertex weight;

[0140] vii. Randomly select a vertex as a seed point according to the established distribution, add it to the seed point list, assign it a weight of 0, mark it as visited and add it to the queue Q;

[0141] viii. Mark all non-seed point vertices as unvisited;

[0142] ix. Repeat steps ii–viii until K seed points are selected.

[0143] Furthermore, the clustering based on diffuse geodesic distance specifically includes:

[0144] The clustering method based on diffuse geodesic distance includes two steps: allocation and seed point update;

[0145] The assignment step assigns vertices to the cluster of the nearest seed point based on geodesic distance:

[0146] i. Initialize the diffuse weighted undirected graph G D : Assign all non-seed point vertex weights to floating-point maximum values ​​and mark them as unassigned;

[0147] ii. Assign the weight of all seed points to 0, mark them as allocated, and add them to the priority queue Q;

[0148] iii. Pop the head vertex v of the queue Q and visit the unvisited adjacent vertices u of v in turn;

[0149] iv. If w v +w e (v,u) <w u , then assign u to the cluster where v is located and mark it as assigned, and update the weight of u to w v +w e (v,u) and join the tail of queue Q, where w v and w u are the weights of vertices v and u, respectively, and w e (v,u) is the weight of the edge connecting vertices v and u;

[0150] v. Iterate steps iii and iv until the priority queue Q is empty;

[0151] Seed point update step The dynamic programming algorithm selects the vertex with the smallest sum of geodesic distances from other vertices in the cluster as the new seed point. For any cluster:

[0152] i. Calculate and record the sum of the vertex weights in the cluster, that is, the sum of the geodesic distances from the seed point s to other vertices in the cluster;

[0153] ii. For the vertex n in the cluster adjacent to the seed point, if the sum of its geodesic distances from the vertices in the cluster is not calculated, execute steps iii–viii;

[0154] iii. Initialize all vertex weights in the cluster to floating point maximum values ​​and mark them as unvisited;

[0155] iv. Assign the weight of vertex n to 0, mark it as visited and add it to the priority queue Q;

[0156] v. Pop the head vertex v of the queue Q, and visit the unvisited adjacent vertices u of v in the cluster in turn;

[0157] vi. If w v +w e (v,u) <w u , then update the weight of u to w v +w e (v,u) <w u And join the tail of queue Q, where w v and w u are the weights of vertices v and u, respectively, and w e (v,u) is the weight of the edge connecting vertices v and u;

[0158] vii. Iterate steps iv–vii until the priority queue Q is empty;

[0159] viii. Sum the weights of the vertices in the cluster and record it as the sum of the geodesics of vertex n from other vertices in the cluster;

[0160] ix. Take the vertex n0 adjacent to the seed point s with the smallest sum of geodesics. If the sum of its geodesics is less than the seed point s, take n0 as the new seed point and execute step ii. Otherwise, terminate the loop.

[0161] Iterate the assignment step and the seed point update step until the seed point does not change or loops in a certain vertex set;

[0162] The voxels corresponding to the vertices in the same cluster form a diffuse supervoxel.

[0163] S103: Taking the diffuse supervoxel as a vertex, generating a weighted undirected graph of supervoxel nerve fiber connections, extracting the maximum cluster of the weighted undirected graph of supervoxel nerve fiber connections, and updating the weighted undirected graph of supervoxel nerve fiber connections according to the maximum cluster to achieve fiber bundle clustering;

[0164] The fiber bundle clustering based on diffuse supervoxels specifically includes:

[0165] The weighted undirected graph G of supervoxel nerve fiber connections is constructed with diffuse supervoxels as vertices. S , if and only if there exists a supervoxel S i and S j When the nerve fiber bundle is connected, there is an edge e(S i ,S j );

[0166] Assume that after S i and S j The fiber bundle sets are T i and T j , the volumes are V i and V j , then the edge e(S i ,S j ) is:

[0167]

[0168] For Figure G S The adjacency matrix of is binarized, and the Otsu algorithm can be used to find the threshold, or the sparsity of the adjacency matrix after binarization |E| / K(K-1) satisfies the given conditions, where |E| is the number of edges after binarization, and K is the number of diffuse supervoxels;

[0169] Use Bron-Kerbosch to extract all the largest groups of the binarized graph, that is, all the largest connected components;

[0170] Take the largest group C with the largest internal connection (i.e. the largest sum of the weights of the internal edges) max ;

[0171] Take the least number of passes through C max The fiber bundles with diffuse supervoxels within ξ% are regarded as clusters, which are deleted from all fiber bundles and the graph G is updated. S The edge weight, preferably, ξ = 60;

[0172] Iterate the above steps until the newly obtained C max The internal connection is less than a given threshold, preferably, the threshold is K×G S Adjacency matrix element mean, K=3.

[0173] S104: registering the diffusion supervoxel to the individual space, segmenting the individual fiber bundles according to the maximum cluster of the individual diffusion supervoxel, performing curve fitting according to the fiber density image of the individual fiber bundle to obtain the fiber bundle center curve, extracting the characteristics of the fiber bundle center curve and performing statistical analysis to achieve segmentation analysis of the fiber bundle.

[0174] The preprocessed dMRI was registered to the constructed dMRI brain template based on the similarity of the water molecule diffusion model (log-Euclidean distance, L2 norm distance of ODF and Jensen-Shannon divergence are optional);

[0175] The diffusion supervoxels in the dMRI brain template space were transformed to the individual space using the inverse deformation of the registration;

[0176] According to the order of the largest cluster of diffuse supervoxels extracted in fiber bundle clustering, when a fiber bundle in the individual space passes through its ξ% diffuse supervoxel, it is considered to belong to the cluster. Preferably, ξ=60.

[0177] The process of obtaining the fiber bundle center curve by curve fitting based on the fiber density image of the individual fiber bundle:

[0178] Up-sample the dMRI brain template and construct fiber bundle density maps of fiber bundle clusters;

[0179] The central curve of the fiber bundle density map is fitted using a non-uniform rational B-spline curve, and the diffusion characteristics of the fiber bundle are mapped to the central curve to obtain the central curve of the fiber bundle.

[0180] Calculation process of fiber bundle center curve:

[0181] Upsampling the dMRI brain template space, preferably, the upsampling factor is 4;

[0182] For any cluster of nerve fiber bundles, the nerve fiber density map was calculated in the upsampled dMRI brain template space, which is the ratio of the number of fiber curves passing through each voxel to the cluster size;

[0183] Use Non-Uniform Rational B-Spline (NURBS) to fit the nerve fiber density map:

[0184]

[0185] Where P i is the control point, N i,p (u) is the B-spline function, Curve(0) and Curve(1) are the center points of the diffuse supervoxels at the two ends of the maximum cluster corresponding to the fiber bundle;

[0186] Let Q0, Q1, …, Qm are the voxel points through which the fiber bundles pass, and their densities are d0, d1, …, d m , NURBS curve fitting uses the least squares method to minimize the objective function:

[0187]

[0188] The NURBS curve is deformed to the individual space using the registered quasi-deformation, and then the diffusion characteristics such as the weighted average FA and GFA of the fiber bundle are mapped to the NURBS curve;

[0189] For the points on the curve, the fiber density map is divided into a Voronoi diagram based on the Euclidean distance, and the diffuse features in the Voronoi diagram are mapped to the corresponding points:

[0190]

[0191] where v k Represents the voxel Q in the Voronoi diagram k Diffusion characteristics. For the NURBS center curve mapping results of the same fiber bundle cluster, the permutation test was used to compare the differences between the groups;

[0192] In the permutation test, for each permutation, the vertices of the curve with significant statistical differences are grouped into clusters according to their adjacency relationships, and the size of the largest cluster is recorded;

[0193] Multiple comparison correction was performed. In the original grouping of the permutation test, if the size of the cluster composed of the vertices of the curve with significant statistical differences was greater than the size of the largest cluster recorded by 95%, the cluster was considered to have statistical differences.

[0194] For clusters with statistical differences, the mean of their diffuse features is taken, and the random forest method is used for feature screening. Then, machine learning classification models such as logistic regression, neural network and support vector machine are constructed.

[0195] Specifically, the present invention is further described below by way of embodiments:

[0196] Embodiment 2: In the second aspect, as Figure 5 As shown, in order to achieve the above purpose, the present invention discloses a neural fiber bundle clustering segmentation system based on diffuse supervoxels, comprising:

[0197] The image acquisition and preprocessing module 11 is used to receive three-dimensional brain magnetic resonance image data, perform preprocessing, obtain processed image data, and register the processed image data using an iterative registration method to obtain a brain template;

[0198] The diffusion supervoxel segmentation module 12 is used to input the brain template into a pre-established water molecule diffusion model, output a diffusion weighted undirected graph, generate a diffusion geodesic distance according to the diffusion weighted undirected graph, and obtain a diffusion supervoxel based on the diffusion geodesic distance segmentation;

[0199] The fiber bundle clustering module 13 is used to generate a weighted undirected graph of supervoxel nerve fiber connections with diffuse supervoxels as vertices, extract the maximum cluster of the weighted undirected graph of supervoxel nerve fiber connections, and update the weighted undirected graph of supervoxel nerve fiber connections according to the maximum cluster to achieve fiber bundle clustering;

[0200] The fiber bundle center curve mapping module 14 is used to align the diffusion supervoxel to the individual space, obtain the individual fiber bundle according to the maximum cluster segmentation of the individual diffusion supervoxel, perform curve fitting based on the fiber density image of the individual fiber bundle to obtain the fiber bundle center curve, extract the characteristics of the fiber bundle center curve and perform statistical analysis to realize the segmentation analysis of the fiber bundle.

[0201] Based on the same inventive concept, the present invention also provides a computer device, which includes: one or more processors, and a memory for storing one or more computer programs; the program includes program instructions, and the processor is used to execute the program instructions stored in the memory. The processor may be a central processing unit (CPU), or other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA) or other programmable logic devices, discrete gates or transistor logic devices, discrete hardware components, etc. It is the computing core and control core of the terminal, which is used to implement one or more instructions, specifically for loading and executing one or more instructions in a computer storage medium to implement the above method.

[0202] It needs to be further explained that, based on the same inventive concept, the present invention also provides a computer storage medium, on which a computer program is stored, and the computer program is executed by a processor to execute the above method. The storage medium can adopt any combination of one or more computer-readable media. The computer-readable medium can be a computer-readable signal medium or a computer-readable storage medium. The computer-readable storage medium can be, for example, but not limited to, an electrical, magnetic, optical, electrical, magnetic, infrared, or semiconductor system, device or device, or any combination of the above. More specific examples (non-exhaustive list) of computer-readable storage media include: an electrical connection with one or more wires, a portable computer disk, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or flash memory), an optical fiber, a portable compact disk read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the above. In the present invention, a computer-readable storage medium can be any tangible medium containing or storing a program, which can be used by an instruction execution system, device or device or used in combination with it.

[0203] In the description of this specification, the description with reference to the terms "one embodiment", "example", "specific example", etc. means that the specific features, structures, materials or characteristics described in conjunction with the embodiment or example are included in at least one embodiment or example of the present disclosure. In this specification, the schematic representation of the above terms does not necessarily refer to the same embodiment or example. Moreover, the specific features, structures, materials or characteristics described can be combined in any one or more embodiments or examples in a suitable manner.

[0204] The above shows and describes the basic principles, main features and advantages of the present disclosure. Those skilled in the art should understand that the present disclosure is not limited by the above embodiments, and the above embodiments and descriptions are only for explaining the principles of the present disclosure. Without departing from the spirit and scope of the present disclosure, the present disclosure may have various changes and improvements, and these changes and improvements fall within the scope of the present disclosure to be protected.

Claims

1. A neural fiber bundle clustering segmentation method based on diffuse supervoxels, characterized in that: The method comprises the following steps: receiving three-dimensional brain magnetic resonance imaging data, performing preprocessing to obtain processed image data, and registering the processed image data using an iterative registration method to obtain a brain template; The brain template is input into a pre-established water molecule diffusion model, and a diffusion weighted undirected graph is outputted. The diffusion geodesic distance is generated according to the diffusion weighted undirected graph, and the diffusion supervoxel is obtained based on the diffusion geodesic distance segmentation; Taking diffuse supervoxels as vertices, a weighted undirected graph of supervoxel nerve fiber connections is generated, the maximum cluster of the weighted undirected graph of supervoxel nerve fiber connections is extracted, and the weighted undirected graph of supervoxel nerve fiber connections is updated according to the maximum cluster to achieve fiber bundle clustering; The diffusion supervoxel is registered to the individual space, and the individual fiber bundles are obtained according to the maximum cluster segmentation of the individual diffusion supervoxel. The fiber density image of the individual fiber bundle is used for curve fitting to obtain the fiber bundle center curve. The characteristics of the fiber bundle center curve are extracted and statistically analyzed to achieve segmentation and analysis of the fiber bundle.

2. The method for clustering and segmenting nerve fiber bundles based on diffuse supervoxels according to claim 1, characterized in that: The process of using the iterative registration method to register the processed image data: Given the processed image data (J1, J2, ..., J n ), the deformations registered to the brain template φ are θ1, θ2, …, θ n , the brain template θ is established by the group registration method: The individual images J1, J2, ..., J n Register to current brain template The deformations are t represents the tth iteration; Pick The non-rigid deformation part of the inverse transformation is averaged as the residual deformation R (t) ; Use residual deformation R (t) Correction deformation in is the corrected deformation; Use the corrected deformation Image J i Transform and take the mean of the transformed image as the brain template for the next iteration n is the number of images; Repeat the above until the difference between two consecutive generated brain templates is less than the preset value ε: The brain template construction process adopts a coarse-to-fine approach. First, the brain template is established using affine transformation, and then nonlinear deformation is used to successively establish low-resolution brain templates and high-resolution brain templates. Use the brain template to build the process, using the ICBM-152 standard brain template as the initial template Establish T1WI brain template; Then use the T1WI brain template to construct the final step T correction deformation The aligned dMRI images are deformed into the T1WI brain template space, and the mean is taken as the initial brain template of dMRI; Based on the initial brain template of dMRI, the dMRI brain template is established by applying the above brain template construction process.

3. The method for clustering and segmenting nerve fiber bundles based on diffuse supervoxels according to claim 1, characterized in that: The process of inputting the brain template into a pre-established water molecule diffusion model and outputting a diffusion weighted undirected graph: Taking white matter voxels as the vertices of the graph, when two voxels have a 26-adjacent relationship, there is an edge e between the corresponding vertices v1 and v2, and the weight of edge e is defined as: W e =d diffusion +γd chamfer where d diffusion Represents the distance between two voxel points of the water molecule diffusion model, d chamfer represents the spatial chamfer distance, γ is the coefficient for adjusting the importance of the spatial chamfer distance; For the diffusion tensor model of water molecule diffusion, the difference is defined as the logarithmic Euclidean distance of the diffusion tensor. Assuming that the diffusion tensors corresponding to v1 and v2 are D1 and D2 respectively, the logarithmic Euclidean distance is: d tensor =tr((log D1-log D2) T (log D1-log D2)) Where tr represents the trace of the matrix; For the directional distribution function model of water molecule diffusion, the ODF model f can represent a linear combination of spherical harmonic functions: where θ and are the polar angle and azimuth of the spherical coordinate system, Y lm is the spherical harmonic function, f lm are the spherical harmonic coefficients, l and m are the degree and order of the spherical harmonic function respectively; Assume that the ODF models of vertices v1 and v2 are f and g respectively, then according to the Bregman divergence, the square of their L2 norm is: Where S represents the unit sphere, Ω represents the spherical angle and θ and are the polar angle and azimuth angle of the spherical coordinate system respectively; since the spherical harmonic function is a set of orthogonal bases, the L2 norm between ODFs is simplified to: In addition, according to the Bregman divergence, the Kullback-Leibler divergence between ODFs is defined as: Then its symmetric non-negative Jensen-Shannon divergence is d KL (f,g) and d KL The mean of (g,f) is: Assume lnf lm and lng lm are the spherical harmonic coefficients of lnf and lng respectively. According to the orthogonality of spherical harmonic functions, the Jensen-Shannon divergence is simplified to: The coefficient γ is used to balance the water molecule diffusion model distance and the space chamfer distance, and is set to γ1 multiplied by the median of the ratio of the water molecule diffusion model distance to the space chamfer distance of the adjacent voxel points.

4. The method for clustering and segmenting nerve fiber bundles based on diffuse supervoxels according to claim 1, characterized in that: The process of generating diffuse geodesic distance according to diffuse weighted undirected graph is as follows: Initialization of seed points based on white matter regions is based on diffusion weighted undirected graph G D Select the seed point; Based on the seed point initialization of the white matter skeleton, the various fractional anisotropies or generalized fractional anisotropy of the white matter are refined by graphics to obtain the white matter skeleton, and then the G corresponding to the white matter skeleton is extracted. D The subgraph G D 'Select a seed point; Assignment step: Calculate the diffuse weighted undirected graph G D The geodesic distance between the vertex and the seed point, assigning the vertex to the cluster corresponding to the seed point with the minimum geodesic distance; Seed point update step: Use the dynamic programming algorithm to select the vertex with the smallest sum of geodesic distances to other vertices in the cluster as the new seed point; Iterate the assignment step and the seed point update step until the seed point does not change or cycles through a few vertex sets; The voxels corresponding to the vertices in the same cluster form a diffuse supervoxel.

5. The method for clustering and segmenting nerve fiber bundles based on diffuse supervoxels according to claim 1, characterized in that: The fiber bundle clustering is performed based on diffuse supervoxels, and the process is as follows: Construct a weighted undirected graph G of supervoxel nerve fiber connections with diffuse supervoxels as vertices S , if and only if there exists a supervoxel S i and S j When the nerve fiber bundle is connected, there is an edge e(S i ,S j ); Assume that after S i and S j The fiber bundle sets are T i and T j , the volumes are V i and V j , then the edge e(S i ,S j ) is: For Figure G S Binarize the adjacency matrix, select the Otsu algorithm to find the threshold, or the sparsity of the adjacency matrix after binarization |E| / K(K-1) satisfies the given conditions, where |E| is the number of edges after binarization and K is the number of diffuse supervoxels; Use Bron-Kerbosch to extract all the largest groups of the binarized graph, which are all the largest connected components; Take the largest clique C with the largest internal connection max ; Take the least number of passes through C max The fiber bundles with diffuse supervoxels within ξ% are regarded as clusters, which are deleted from all fiber bundles and the graph G is updated. S The weight of the edge, ξ is the preset ratio; Iterate the above steps until the newly obtained C max The internal connection is less than a given threshold, where the threshold is K×G S Adjacency matrix element mean.

6. The method for clustering and segmenting nerve fiber bundles based on diffuse supervoxels according to claim 1, characterized in that: The process of performing curve fitting based on the fiber density image of the individual fiber bundle to obtain the fiber bundle center curve is as follows: Up-sample the dMRI brain template and construct fiber bundle density maps of fiber bundle clusters; The central curve of the fiber bundle density map is fitted using a non-uniform rational B-spline curve, and the diffusion characteristics of the fiber bundle are mapped to the central curve to obtain the central curve of the fiber bundle.

7. The method for clustering and segmenting nerve fiber bundles based on diffuse supervoxels according to claim 6, characterized in that: The calculation process of the fiber bundle center curve is as follows: For any cluster of nerve fiber bundles, the nerve fiber density map was calculated in the upsampled dMRI brain template space, which was the ratio of the number of fiber curves passing through each voxel to the cluster size; Use non-uniform rational B-spline curve to fit the nerve fiber density map: Where P i is the control point, N i,p (u) is the B-spline function, u takes values ​​between 0 and 1, Curve(0) and Curve(1) are the center points of the diffuse supervoxels at the two ends of the maximum cluster corresponding to the fiber bundle; Let Q0, Q1, …, Q m are the voxel points through which the fiber bundles pass, and their densities are d0, d1, …, d m , NURBS curve fitting uses the least squares method to minimize the objective function: The NURBS curve is deformed to the individual space using the registered quasi-deformation, and then the diffusion characteristics such as the weighted average FA and GFA of the fiber bundle are mapped to the NURBS curve; For the points on the curve, the fiber density map is divided into a Voronoi diagram based on the Euclidean distance, and the diffuse features in the Voronoi diagram are mapped to the corresponding points: Σd k v k / Σd k where v k Represents the voxel Q in the Voronoi diagram k diffusion characteristics.

8. A diffuse supervoxel-based neural fiber bundle clustering segmentation system, which adopts the diffuse supervoxel-based neural fiber bundle clustering segmentation method according to any one of claims 1 to 7, characterized in that: include: An image acquisition and preprocessing module is used to receive three-dimensional brain magnetic resonance image data, perform preprocessing, obtain processed image data, and register the processed image data using an iterative registration method to obtain a brain template; A diffusion supervoxel segmentation module is used to input the brain template into a pre-established water molecule diffusion model, output a diffusion weighted undirected graph, generate a diffusion geodesic distance based on the diffusion weighted undirected graph, and obtain a diffusion supervoxel based on the diffusion geodesic distance segmentation; A fiber bundle clustering module is used to generate a weighted undirected graph of supervoxel nerve fiber connections with diffuse supervoxels as vertices, extract the maximum cluster of the weighted undirected graph of supervoxel nerve fiber connections, and update the weighted undirected graph of supervoxel nerve fiber connections according to the maximum cluster to achieve fiber bundle clustering; The fiber bundle center curve mapping module is used to align the diffusion supervoxel to the individual space, obtain the individual fiber bundle according to the maximum cluster segmentation of the individual diffusion supervoxel, perform curve fitting based on the fiber density image of the individual fiber bundle to obtain the fiber bundle center curve, extract the characteristics of the fiber bundle center curve and perform statistical analysis to realize the segmentation analysis of the fiber bundle.

9. A terminal device, comprising a memory, a processor, and a computer program stored in the memory and capable of running on the processor, characterized in that: The memory stores a computer program that can be run on the processor. When the processor loads and executes the computer program, the neural fiber bundle clustering segmentation method based on diffuse supervoxels described in any one of claims 1 to 7 is adopted.

10. A computer-readable storage medium having a computer program stored therein, characterized in that: When the computer program is loaded and executed by the processor, the neural fiber bundle clustering segmentation method based on diffuse supervoxels according to any one of claims 1 to 7 is adopted.

Images