Diffusion super voxel based neural fiber bundle clustering segmentation method and system

Abstract

This invention discloses a method and system for neural fiber tract clustering and segmentation based on diffuse hypervoxels, relating to the fields of medical image analysis and artificial intelligence. The method includes: receiving and preprocessing three-dimensional brain magnetic resonance imaging (MRI) data; establishing a population-matched personalized MRI brain template using an iterative registration method; constructing a diffusion-weighted undirected graph and segmenting diffuse hypervoxels based on diffusion geodesic distance; clustering fiber tracts based on the maximum clusters of the neural fiber connectivity graph of the diffuse hypervoxels; registering the diffuse hypervoxels to an individual space and using the maximum clusters to achieve personalized fiber tract segmentation; constructing a center curve for the fiber tract clusters and mapping diffusion features to this curve; establishing a statistical analysis method for the center curve based on permutation tests and building a machine learning classification model.

Description

Technical Field

[0001] This invention relates to the fields of medical image analysis and artificial intelligence technology, specifically to a method and system for clustering and segmenting neural fiber bundles based on diffuse hypervoxels. Background Technology

[0002] Diffusion magnetic resonance imaging (dMRI) uses a gradient field to alter 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 spectroscopic imaging (DSI), and high angular resolution diffusion imaging (HARDI). Fiber tractor methods establish diffusion tensor models or orientation distribution functions (ODFs) of water molecule diffusion from dMRI, and then reconstruct fiber curves describing the course and orientation of white matter fiber bundles.

[0003] Fiber tractography results can reflect the damage to white matter fiber tracts caused by disease, thereby improving diagnostic accuracy. Since multiple fiber tracts connect different brain regions, it is necessary to cluster the fiber tractography results to construct a white matter fiber tractograph, and then segment individual white matter fiber tracts based on the fiber tractograph. Current fiber tract clustering methods are mostly based on the distance between fiber curves, such as Hausdorff distance and mean shortest distance, which have high computational complexity, require downsampling of fiber tractography results, and are not suitable for constructing fiber tractographs based on large-sample population dMRI image data, thus limiting the development of related fields. Summary of the Invention

[0004] To address the shortcomings mentioned in the background art, the present invention aims to provide a method and system for clustering and segmenting neural fiber bundles based on diffuse hypervoxels, which transforms the white matter fiber bundle clustering problem into a diffuse hypervoxel clustering problem, thereby reducing computational complexity and increasing computational speed.

[0005] Firstly, the objective of this invention can be achieved through the following technical solution: a neural fiber bundle clustering and segmentation method based on diffuse hypervoxels, the method comprising the following steps:

[0006] The system receives three-dimensional brain magnetic resonance imaging data, performs preprocessing to obtain processed image data, and uses an iterative registration method to register the processed image data 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 output. The diffusion geodesic distance is generated based on the diffusion weighted undirected graph, and the diffusion hypervoxel is obtained by segmentation based on the diffusion geodesic distance.

[0008] Using diffuse hypervoxels as vertices, a weighted undirected graph of hypervoxel neural fiber connections is generated. The largest clique in the weighted undirected graph of hypervoxel neural fiber connections is extracted, and the weighted undirected graph of hypervoxel neural fiber connections is updated based on the largest clique to achieve fiber bundle clustering.

[0009] Diffused hypervoxels are registered to individual space, and individual fiber bundles are obtained by segmenting based on the maximum cluster of individual diffuse hypervoxels. Curve fitting is performed on the fiber density image of the individual fiber bundle to obtain the fiber bundle center curve. The features of the fiber bundle center curve are extracted and statistically analyzed to achieve fiber bundle segmentation analysis.

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

[0011] Given the processed image data (J1, J2, ..., J... n The deformations of the brain template φ registered to it are θ1, θ2, ..., θ n The brain template θ is established using a group registration method:

[0012] Individual images J1, J2, ..., J n Register to the current brain template The deformations are respectively t represents the t-th iteration;

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

[0014] Using residual deformation R (t) Correction deformation in The deformation is corrected;

[0015] Using the corrected deformation For image J i The transformation is performed, and the mean of the transformed image is used 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, affine transformation is used to build the brain template, and then nonlinear deformation is used to build low-resolution brain templates and high-resolution brain templates in sequence.

[0019] The brain template construction process was used, with the ICBM-152 standard brain template as the initial template. Establish T1WI brain template;

[0020] Then, the final step in the T1WI brain template construction process was to correct the deformation. The aligned dMRI images were deformed into the T1WI brain template space, and the mean value was taken as the initial brain template for dMRI.

[0021] The initial brain template is based on dMRI, and the above brain template construction process is applied to build a dMRI brain template.

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

[0023] Using white matter voxels as vertices in a graph, when two voxels are 26 adjacent, there is an edge e between the corresponding vertices v1 and v2. The weight of edge e is defined as follows:

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

[0025] Where d diffusion d represents the distance in the two-voxel point water molecule diffusion model. chamfer Indicates the spatial chamfer distance;

[0026] For the dispersion tensor model of water molecule diffusion, the difference is defined as the logarithmic Euclidean distance of the dispersion tensors. Let the dispersion tensors corresponding to v1 and v2 be D1 and D2, respectively, then 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 be represented as a linear combination of spherical harmonic functions:

[0030]

[0031] Where θ and These are the polar angle and azimuth angle of the spherical coordinate system, respectively. lm f is a spherical harmonic function. lm Let be the spherical harmonic coefficients, and l and m be the degree and order of the spherical harmonic function, respectively;

[0032] Let f and g be the ODF models of v1 and v2, respectively. Then, according to the Bregman divergence, the squares of their L2 norms are:

[0033]

[0034] Where S represents the unit sphere, Ω represents the spherical angle, and θ and Let be the polar angle and azimuth angle of the spherical coordinate system, respectively; since the spherical harmonics are an orthogonal basis, the L2 norm between ODFs simplifies to:

[0035]

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

[0037]

[0038] Then its symmetric nonnegative Jensen-Shannon divergence is d KL (f,g) and d KL Mean of (g,f):

[0039]

[0040] Let lnf lm and lng lm Let lnf and lng be the spherical harmonic coefficients, respectively. Then, based on the orthogonality of spherical harmonic functions, the Jensen-Shannon divergence simplifies to:

[0041]

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

[0043] In conjunction with the first aspect, in some implementations of the first aspect, the method further includes: generating diffuse geodesic distances based on the diffuse weighted undirected graph, and initializing seed points using the KMeans++ algorithm based on the geodesic distances between vertices in the diffuse weighted undirected graph, as follows:

[0044] Seed point initialization based on white matter region is based on diffusion-weighted undirected graph G. D Select a seed point;

[0045] The white matter skeleton is obtained by initializing seed points based on the white matter skeleton and performing graphical refinement operations on the fractional anisotropy or generalized anisotropy of the white matter. Then, the G corresponding to the white matter skeleton is extracted. D subgraph G′ D Select a seed point;

[0046] Allocation steps: Calculate the diffuse weighted undirected graph G D The geodesic distance between a vertex and a seed point is used to assign the vertex to the cluster corresponding to the seed point with the smallest geodesic distance.

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

[0048] Iterative allocation and seed point update steps continue until the seed point remains unchanged or loops through a few vertex sets; voxels corresponding to vertices within the same cluster form a diffuse supervoxel.

[0049] In conjunction with the first aspect, in some implementations of the first aspect, the method further includes: the fiber bundle clustering is based on diffuse hypervoxels, as follows:

[0050] A weighted undirected graph G is constructed using diffuse supervoxels as vertices to connect supervoxel neural fibers. S If and only if there exists a supervoxel S that has passed through it simultaneously i and S j When nerve fiber bundles are present, there is a boundary e(S) between them. i ,S j );

[0051] Let S pass through i and S j The fiber bundle sets are T i and T j The volumes are V i and V j Then edge e(S) i ,S j The weights of ) are:

[0052]

[0053] For graph G SThe adjacency matrix is ​​binarized. The threshold can be found using the Otsu algorithm, or the sparsity of the adjacency matrix after binarization, |E| / K(K-1), can satisfy a given condition, where |E| is the number of edges after binarization and K is the number of diffuse hypervoxels.

[0054] Use Bron-Kerbosch to extract all the maximum cliques in the binarized graph, which are then used to extract all the maximum connected components.

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

[0056] Take the minimum number of passes through C max The fiber bundles of diffuse supervoxels within the ξ% region are treated as clusters, removed from all fiber bundles, and graph G is updated. S The weight of the edge, ξ, is a preset ratio;

[0057] Iterate the above steps until a new C is obtained. max The number of internal connections is less than a given threshold, where the threshold is K×G. S Mean of elements in the adjacency matrix.

[0058] In conjunction with the first aspect, in some implementations of the first aspect, the method further includes: the process of obtaining the fiber bundle center curve by curve fitting based on the fiber density image of the individual fiber bundle:

[0059] Upsample the dMRI brain template and construct a fiber tract density map of the fiber tract clusters;

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

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

[0062] For any cluster of nerve fibers, the nerve fiber density map is 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.

[0063] Nerve fiber density maps were fitted using non-uniform rational B-spline curves.

[0064]

[0065] Where P i N is the control point. i,p (u) is a B-spline function, with u taking values ​​between 0 and 1. Curve(0) and Curve(1) are the center points of the diffuse hypervoxels at both ends of the maximum cluster corresponding to the fiber bundle.

[0066] Let Q0, Q1, ..., Q m The voxel points through which the fiber bundle passes have densities d0, d1, ..., d m NURBS curve fitting uses the least squares method to minimize the objective function:

[0067]

[0068] The NURBS curves are deformed to the individual space using the registration-based pseudo-deformation, and then the diffuse features such as the weighted average FA and GFA of the fiber bundles are mapped to the NURBS curves.

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

[0070]

[0071] Where v k Voxel Q in Voronoi diagram k The diffusion characteristics.

[0072] Secondly, in order to achieve the above objectives, the present invention discloses a neural fiber bundle clustering and segmentation system based on diffuse hypervoxels, comprising:

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

[0074] The diffusion hypervoxel segmentation module is used to input the brain template into a pre-established water molecule diffusion-based model and output a diffusion-weighted undirected graph. The diffusion geodesic distance is generated based on the diffusion-weighted undirected graph, and the diffusion hypervoxel is segmented based on the diffusion geodesic distance.

[0075] The fiber bundle clustering module is used to generate a weighted undirected graph of supervoxel neural fiber connections with diffuse supervoxels as vertices, extract the maximum clique of the weighted undirected graph of supervoxel neural fiber connections, and update the weighted undirected graph of supervoxel neural fiber connections based on the maximum clique to achieve fiber bundle clustering.

[0076] The fiber bundle center curve mapping module is used to register diffuse hypervoxels to individual space, segment individual fiber bundles based on the maximum cluster of individual diffuse hypervoxels, perform curve fitting based on the fiber density image of individual fiber bundles to obtain fiber bundle center curves, extract the features of fiber bundle center curves and perform statistical analysis to realize fiber bundle segmentation analysis.

[0077] In another aspect of the present invention, in order to achieve the above-mentioned objective, 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. The memory stores the computer program capable of running on the processor. When the processor loads and executes the computer program, it employs the neural fiber bundle clustering and segmentation method based on diffuse hypervoxels as described above.

[0078] In another aspect of the present invention, in order to achieve the above-mentioned objective, a computer-readable storage medium is disclosed, wherein the computer-readable storage medium stores a computer program, characterized in that, when the computer program is loaded and executed by a processor, it employs the neural fiber bundle clustering segmentation method based on diffuse hypervoxels as described above.

[0079] The beneficial effects of this invention are:

[0080] This invention transforms the complex fiber bundle clustering problem into an easily solvable diffuse supervoxel clustering problem, which has the advantages of fast computation speed and low memory requirements. It can handle large-scale fiber curves and provides a feasible solution for constructing fiber bundle maps of large sample populations. It segments fiber bundles according to the spatial relationship of diffuse supervoxels, which has the advantages of fast computation speed and high robustness, and provides a new solution for statistical analysis and classification model training of specific fiber bundles. Attached Figure Description

[0081] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, for those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0082] Figure 1 This is a flowchart of the method of the present invention;

[0083] Figure 2 This is a flowchart of the dMRI brain template construction in an embodiment of the present invention;

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

[0085] Figure 4 This is a flowchart of fiber bundle clustering based on diffuse hypervoxels in an embodiment of the present invention;

[0086] Figure 5 This is a schematic diagram of the system structure of the present invention. Detailed Implementation

[0087] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0088] Example 1:

[0089] The following is a description of the relevant terms used in the embodiments of this application:

[0090] Nerve fiber bundle: In white matter, a bundle of nerve fibers that have similar origin, termination, course and function is formed by grouping together.

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

[0092] S101: Receives three-dimensional brain magnetic resonance imaging data, performs preprocessing to obtain processed image data, and uses an iterative registration method to register the processed image data to obtain a brain template;

[0093] Preprocessing steps include removing non-brain tissue, gray and white matter segmentation, registration, water molecule diffusion model reconstruction, and fiber bundle tracing.

[0094] The process of registering 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 of the brain template φ registered to it are θ1, θ2, ..., θ n The brain template θ can be established using group registration methods:

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

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

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

[0099] iv. Use the modified deformation For image J iThe transformation is performed, and the mean of the transformed image is used as the brain template for the next iteration.

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

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

[0102] Using the brain template construction process described above, with the ICBM-152 standard brain template as the initial template. Establish T1WI brain template;

[0103] Then, the final step in the T1WI brain template construction process was to correct the deformation. The aligned dMRI images were deformed into the T1WI brain template space, and the mean value was taken as the initial brain template for dMRI.

[0104] Based on the initial brain template from dMRI, the above-described brain template construction process was used to build a dMRI brain template.

[0105] S102: Input the brain template into the pre-established water molecule diffusion model and output a diffusion weighted undirected graph. Generate diffusion geodesic distance based on the diffusion weighted undirected graph and obtain diffusion hypervoxels based on the diffusion geodesic distance segmentation.

[0106] like Figure 3 As shown, diffusion-weighted undirected graph G is constructed based on the difference between the spatial adjacency relationship of voxel points in the white matter brain region and the water molecule diffusion model. D Then, seed points (i.e., the central voxels of clusters) are initialized in the white matter region or white matter skeleton region using the KMeans++ method based on diffusion distance. Finally, the white matter region is over-segmented into diffusion super voxels using a clustering method based on diffusion geodesic distance.

[0107] Diffuse weighted undirected graph G D The construction, specifically includes:

[0108] Using white matter voxels as vertices in the graph, when two voxels are 26-adjacent (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 follows:

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

[0110] Where d diffusiond represents the distance in the two-voxel point water molecule diffusion model. chamfer Indicates the spatial chamfer distance;

[0111] For the dispersion tensor model of water molecule diffusion, the difference is defined as the logarithmic Euclidean distance of the dispersion tensors. Let the dispersion tensors corresponding to v1 and v2 be D1 and D2, respectively, then 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 be represented as a linear combination of spherical harmonic functions:

[0115]

[0116] Where θ and These are the polar angle and azimuth angle of the spherical coordinate system, respectively. lm f is a spherical harmonic function. lm Let be the spherical harmonic coefficients, and l be the order of the spherical harmonic function;

[0117] Let f and g be the ODF models of v1 and v2, respectively. Then, according to the Bregman divergence, the squares of their L2 norms are:

[0118]

[0119] Where S represents the unit sphere, Ω represents the spherical angle, and Since the spherical harmonics form an orthogonal basis, the L2 norm between ODFs simplifies to:

[0120]

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

[0122]

[0123] Then its symmetric nonnegative Jensen-Shannon divergence is d KL (f,g) and d KL Mean of (g,f):

[0124]

[0125] Let lnflm and lng lm Let lnf and lng be the spherical harmonic coefficients, respectively. Then, based on the orthogonality of spherical harmonic functions, the Jensen-Shannon divergence can be simplified as:

[0126]

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

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

[0129] Generate diffuse geodesic distances from the diffuse weighted undirected graph. Initialize seed points using the KMeans++ algorithm based on the geodesic distances between vertices in the diffuse weighted undirected graph. The process is as follows:

[0130] Seed point initialization based on white matter region is based on diffusion-weighted undirected graph G. D Select a seed point;

[0131] Based on seed point initialization of the white matter skeleton, a graphical refinement operation is performed on the fractional anisotropy (FA) or generalized fractional anisotropy (GFA) of the white matter to obtain the white matter skeleton. Then, the corresponding GFA is extracted from the white matter skeleton. D subgraph G′ D Select a seed point.

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

[0133] i. Initialize G, assign all vertex weights to floating-point maxima, 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 from 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 If u is visited, then mark u as visited and update its weight to w. v +w e(v,u) and add it to the tail of queue Q, where w v and w u The weights of vertices v and u are respectively, w e (v,u)

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

[0138] v. Repeat steps iii and iv until queue Q is empty, and the vertex weight is its 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 weights;

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

[0141] viii. Mark all non-seed 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] Clustering methods based on diffuse geodesic distance include two steps: assignment and seed point update.

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

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

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

[0148] iii. Pop the head vertex v from 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 If u is assigned to the cluster where v is located and marked as assigned, then the weight of u is updated to w. v +w e (v,u) and add it to the tail of queue Q, where w v and w u The weights of vertices v and u are respectively, w e (v,u) represents the weight of the edge connecting vertices v and u;

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

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

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

[0153] ii. For a vertex n in the cluster adjacent to the seed point, if the sum of its geodesic distances to other vertices in the cluster has not been calculated, proceed to steps iii–viii.

[0154] iii. Initialize all vertex weights within the cluster to non-floating-point maxima and mark them as unvisited;

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

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

[0157] vi. If w v +w e (v,u) <w u Then the weight of u will be updated to w. v +w e (v,u) <w u And add it to the tail of queue Q, where w v and w u The weights of vertices v and u are respectively, w e (v,u) represents the weight of the edge connecting vertices v and u;

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

[0159] viii. Sum the vertices within the cluster and record them as the sum of the geodesic distances of vertex n from other vertices within the cluster;

[0160] ix. Select the vertex n0 that is adjacent to the seed point s and has the minimum sum of geodesics. If the sum of its geodesics is less than that of the seed point s, select n0 as the new seed point and execute step ii. Otherwise, terminate the loop.

[0161] Iteratively execute the allocation step and the seed point update step until the seed point no longer changes or loops through a few vertex sets;

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

[0163] S103: Using diffuse hypervoxels as vertices, generate a weighted undirected graph of hypervoxel neural fiber connections, extract the largest clique of the weighted undirected graph of hypervoxel neural fiber connections, update the weighted undirected graph of hypervoxel neural fiber connections based on the largest clique, and realize fiber bundle clustering;

[0164] Fiber bundle clustering based on diffuse hypervoxels specifically includes:

[0165] A weighted undirected graph G is constructed using diffuse hypervoxels as vertices to connect hypervoxel neural fibers. S If and only if there exists a supervoxel S that has passed through it simultaneously i and S j When nerve fiber bundles are present, there is a boundary e(S) between them. i ,S j );

[0166] Let S pass through i and S j The fiber bundle sets are T i and T j The volumes are V i and V j Then edge e(S) i ,S j The weights of ) are:

[0167]

[0168] For graph G S The adjacency matrix is ​​binarized. The threshold can be found using the Otsu algorithm, or the sparsity of the adjacency matrix after binarization, |E| / K(K-1), can satisfy a given condition, where |E| is the number of edges after binarization and K is the number of diffuse hypervoxels.

[0169] Use Bron-Kerbosch to extract all the maximum cliques, i.e. all the maximum connected components, of the binarized graph;

[0170] Find the maximum clique C with the maximum internal connectivity (i.e., the maximum sum of the weights of the internal edges). max ;

[0171] Take the minimum number of passes through C max The fiber bundles of diffuse supervoxels within the ξ% region are treated as clusters, removed from all fiber bundles, and graph G is updated. S The preferred edge weight is ξ = 60.

[0172] Iterate the above steps until a new C is obtained. max The number of internal connections is less than a given threshold, preferably K×G. S The mean of the elements in the adjacency matrix, K=3.

[0173] S104: Register the diffuse hypervoxel to the individual space, segment the individual fiber bundle based on the maximum cluster of the individual diffuse hypervoxel, perform curve fitting based on the fiber density image of the individual fiber bundle to obtain the fiber bundle center curve, extract the features of the fiber bundle center curve and perform statistical analysis to realize the 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 (optional log Euclidean distance, L2 norm distance of ODF and Jensen-Shannon divergence).

[0175] The inverse deformation of the registration was used to transform the diffuse hypervoxels in the dMRI brain template space to the individual space;

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

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

[0178] Upsample the dMRI brain template and construct a fiber tract density map of the fiber tract clusters;

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

[0180] The calculation process of the fiber bundle center curve:

[0181] Upsample the dMRI brain template space, preferably by a factor of 4;

[0182] For any cluster of nerve fibers, calculate the nerve fiber density map 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] The nerve fiber density map was fitted using non-uniform rational B-spline (NURBS) curves:

[0184]

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

[0186] Let Q0, Q1, ..., Qm The voxel points through which the fiber bundle passes have densities d0, d1, ..., d m NURBS curve fitting uses the least squares method to minimize the objective function:

[0187]

[0188] The NURBS curves are deformed to the individual space using the registration-based pseudo-deformation, and then the diffuse features such as the weighted average FA and GFA of the fiber bundles are mapped to the NURBS curves.

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

[0190]

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

[0192] In the permutation test, for curve vertices that show statistically significant differences in each permutation, clusters are formed based on their adjacency relationships, and the size of the largest cluster is recorded.

[0193] After performing multiple comparison correction, if the size of the cluster formed by the vertices of the curves with statistically significant differences in the original groupings of the permutation test is greater than the size of the largest cluster recorded in 95% of the records, the cluster is considered to have statistical differences.

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

[0195] Specifically, the present invention will be further illustrated below through embodiments:

[0196] Example 2: Second aspect, such as Figure 5 As shown, in order to achieve the above objectives, this invention discloses a neural fiber tract clustering and segmentation system based on diffuse hypervoxels, comprising:

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

[0198] The diffusion hypervoxel 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 diffusion geodesic distances based on the diffusion weighted undirected graph, and segment diffusion hypervoxels based on the diffusion geodesic distances.

[0199] The fiber bundle clustering module 13 is used to generate a weighted undirected graph of supervoxel neural fiber connections with diffuse supervoxels as vertices, extract the largest clique of the weighted undirected graph of supervoxel neural fiber connections, and update the weighted undirected graph of supervoxel neural fiber connections according to the largest clique to realize fiber bundle clustering.

[0200] The fiber bundle center curve mapping module 14 is used to register diffuse hypervoxels to individual space, segment individual fiber bundles based on the maximum cluster of individual diffuse hypervoxels, perform curve fitting based on the fiber density image of individual fiber bundles to obtain fiber bundle center curves, extract the features of fiber bundle center curves and perform statistical analysis to realize fiber bundle segmentation analysis.

[0201] Based on the same inventive concept, this invention also provides a computer device, comprising: one or more processors, and a memory for storing one or more computer programs; the programs include program instructions, and the processor executes the program instructions stored in the memory. The processor may be a Central Processing Unit (CPU), or other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. It is the computing and control core of the terminal, used to implement one or more instructions, specifically for loading and executing one or more instructions stored in a computer storage medium to implement the above-described method.

[0202] It should be further explained that, based on the same inventive concept, the present invention also provides a computer storage medium storing a computer program, which, when executed by a processor, performs the above-described method. This storage medium can be 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, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any combination thereof. More specific examples of computer-readable storage media (a non-exhaustive list) include: an electrical connection having one or more wires, a portable computer disk, a hard disk, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fiber, portable compact disk read-only memory (CD-ROM), optical storage device, magnetic storage device, or any suitable combination thereof. In the present invention, the computer-readable storage medium can be any tangible medium containing or storing a program that can be used by or in conjunction with an instruction execution system, apparatus, or device.

[0203] In the description of this specification, references to terms such as "an embodiment," "example," "specific example," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of this disclosure. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples.

[0204] The foregoing has shown and described the basic principles, main features, and advantages of this disclosure. Those skilled in the art should understand that this disclosure is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of this disclosure. Various changes and modifications can be made to this disclosure without departing from its spirit and scope, and all such changes and modifications fall within the scope of this disclosure as claimed.

Claims

1. A neural fiber bundle clustering and segmentation method based on diffuse hypervoxels, characterized in that, The method includes the following steps: The system receives three-dimensional brain magnetic resonance imaging data, performs preprocessing to obtain processed image data, and uses an iterative registration method to register the processed image data 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 output. The diffusion geodesic distance is generated based on the diffusion weighted undirected graph, and the diffusion hypervoxel is obtained by segmentation based on the diffusion geodesic distance. Using diffuse hypervoxels as vertices, a weighted undirected graph of hypervoxel neural fiber connections is generated. The largest clique in the weighted undirected graph of hypervoxel neural fiber connections is extracted, and the weighted undirected graph of hypervoxel neural fiber connections is updated based on the largest clique to achieve fiber bundle clustering. Diffused hypervoxels are registered to individual space, and individual fiber bundles are obtained by segmenting based on the maximum cluster of individual diffuse hypervoxels. Curve fitting is performed on the fiber density image of the individual fiber bundle to obtain the fiber bundle center curve. The features of the fiber bundle center curve are extracted and statistically analyzed to achieve fiber bundle segmentation analysis.

2. The neural fiber bundle clustering and segmentation method based on diffuse hypervoxels according to claim 1, characterized in that, The process of registering the processed image data using an iterative registration method: Given the processed image data (J1, J2, ..., J... n The deformations of the brain template φ registered to it are θ1, θ2, ..., θ n The brain template θ is established using a group registration method: Individual images J1, J2, ..., J n Register to the current brain template The deformations are respectively t represents the t-th iteration; Pick The non-rigid deformation portion of the inverse transformation is averaged as the residual deformation R. (t) ; Using residual deformation R (t) Correction deformation in The deformation is corrected; Using the corrected deformation For image J i The transformation is performed, and the mean of the transformed image is used 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, affine transformation is used to build the brain template, and then nonlinear deformation is used to build low-resolution brain templates and high-resolution brain templates in sequence. The brain template construction process was used, with the ICBM-152 standard brain template as the initial template. Establish T1WI brain template; Then, the final step of the T-wave correction deformation process during the T1WI brain template construction was performed. The aligned dMRI images were deformed into the T1WI brain template space, and the mean value was taken as the initial brain template for dMRI. The initial brain template is based on dMRI, and the above brain template construction process is applied to build a dMRI brain template.

3. The neural fiber bundle clustering and segmentation method based on diffuse hypervoxels according to claim 1, characterized in that, The process of inputting a brain template into a pre-established water molecule diffusion model and outputting a diffusion-weighted undirected graph: Using white matter voxels as vertices in a graph, when two voxels are 26 adjacent, there is an edge e between the corresponding vertices v1 and v2. The weight of edge e is defined as follows: W e =d diffusion +γd chamfer Where d diffusion d represents the distance in the two-voxel point water molecule diffusion model. chamfer This represents the spatial chamfer distance, where γ is a coefficient that adjusts the importance of the spatial chamfer distance. For the dispersion tensor model of water molecule diffusion, the difference is defined as the logarithmic Euclidean distance of the dispersion tensors. Let the dispersion tensors corresponding to v1 and v2 be D1 and D2, respectively, then 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 be represented as a linear combination of spherical harmonic functions: Where θ and These are the polar angle and azimuth angle of the spherical coordinate system, respectively. lm f is a spherical harmonic function. lm Let be the spherical harmonic coefficients, and l and m be the degree and order of the spherical harmonic function, respectively; Let f and g be the ODF models of vertices v1 and v2, respectively. Then, according to the Bregman divergence, the squares of their L2 norms are: Where S represents the unit sphere, Ω represents the spherical angle, and θ and Let be the polar angle and azimuth angle of the spherical coordinate system, respectively; since the spherical harmonics are an orthogonal basis, the L2 norm between ODFs simplifies to: Furthermore, according to the Bregman divergence, the Kullback-Leibler divergence between ODFs is defined as: Then its symmetric nonnegative Jensen-Shannon divergence is d KL (f,g) and d KL Mean of (g,f): Let lnf lm and lng lm Let lnf and lng be the spherical harmonic coefficients, respectively. Then, based on the orthogonality of spherical harmonic functions, the Jensen-Shannon divergence simplifies to: The coefficient γ is used to balance the distance of the water molecule diffusion model and the spatial chamfer distance, and is set to γ1 multiplied by the median of the ratio of the distance of the water molecule diffusion model to the spatial chamfer distance between adjacent voxels.

4. The neural fiber bundle clustering and segmentation method based on diffuse hypervoxels according to claim 1, characterized in that, The process of generating diffuse geodesic distances from the diffuse weighted undirected graph, and initializing seed points using the KMeans++ algorithm based on the geodesic distances between vertices in the graph, is as follows: Seed point initialization based on white matter region is based on diffusion-weighted undirected graph G. D Select a seed point; The white matter skeleton is obtained by initializing seed points based on the white matter skeleton and performing graphical refinement operations on the fractional anisotropy or generalized anisotropy of the white matter. Then, the G corresponding to the white matter skeleton is extracted. D subgraph G D Select seed point; Allocation steps: Calculate the diffuse weighted undirected graph G D The geodesic distance between a vertex and a seed point is used to assign the vertex to the cluster corresponding to the seed point with the smallest geodesic distance. Seed point update steps: Use dynamic programming algorithm to select the vertex in the cluster with the smallest sum of geodesic distances to other vertices as the new seed point; Iterate through the allocation steps and seed point update steps until the seed point no longer changes or loops through a few vertex sets; The voxels corresponding to the vertices within the same cluster form a diffuse supervoxel.

5. The neural fiber bundle clustering and segmentation method based on diffuse hypervoxels according to claim 1, characterized in that, The fiber bundle clustering is based on diffuse supravoxels, and the process is as follows: A weighted undirected graph G is constructed using diffuse supervoxels as vertices to connect supervoxel neural fibers. S If and only if there exists a supervoxel S that has passed through it simultaneously i and S j When nerve fiber bundles are present, there is a boundary e(S) between them. i ,S j ); Let S pass through i and S j The fiber bundle sets are T i and T j The volumes are V i and V j Then edge e(S) i ,S j The weights of ) are: For graph G S The adjacency matrix is ​​binarized, and the threshold is found by the Otsu algorithm, or the sparsity of the adjacency matrix after binarization, |E| / K(K-1), satisfies the given condition, where |E| is the number of edges after binarization and K is the number of diffuse hypervoxels. Use Bron-Kerbosch to extract all the maximum cliques in the binarized graph, which are then used to extract all the maximum connected components. Take the largest clique C with the largest internal connectivity. max ; Take the minimum number of passes through C max The fiber bundles of diffuse supervoxels within the ξ% region are treated as clusters, removed from all fiber bundles, and graph G is updated. S The weight of the edge, ξ, is a preset ratio; Iterate the above steps until a new C is obtained. max The number of internal connections is less than a given threshold, where the threshold is K×G. S Mean of elements in the adjacency matrix.

6. The neural fiber bundle clustering and segmentation method based on diffuse hypervoxels according to claim 1, characterized in that, The process of obtaining the fiber bundle center curve by curve fitting based on the fiber density image of an individual fiber bundle: Upsample the dMRI brain template and construct a fiber tract density map of the fiber tract clusters; The fiber bundle density map center curve is fitted using a non-uniform rational B-spline curve, and the diffusion characteristics of the fiber bundle are mapped onto the center curve to obtain the fiber bundle center curve.

7. The neural fiber bundle clustering and segmentation method based on diffuse hypervoxels according to claim 6, characterized in that, The calculation process of the fiber bundle center curve: For any cluster of nerve fibers, the nerve fiber density map is 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. Nerve fiber density maps were fitted using non-uniform rational B-spline curves. Where P i N is the control point. i,p (u) is a B-spline function, with u taking values ​​between 0 and 1. Curve(0) and Curve(1) are the center points of the diffuse hypervoxels at both ends of the maximum cluster corresponding to the fiber bundle. Let Q0, Q1, ..., Q m The voxel points through which the fiber bundle passes have densities d0, d1, ..., d m NURBS curve fitting uses the least squares method to minimize the objective function: The NURBS curves are deformed to the individual space using the registration-based pseudo-deformation, and then the diffuse features such as the weighted average FA and GFA of the fiber bundles are mapped to the NURBS curves. For points on the curve, the fiber density map is divided into Voronoi diagrams based on Euclidean distance, and the diffuse features within the Voronoi diagrams are mapped to the corresponding points: Σd k v k / Σd k Where v k Voxel Q in Voronoi diagram k The diffusion characteristics.

8. A neural fiber tract clustering and segmentation system based on diffuse hypervoxels, employing the neural fiber tract clustering and segmentation method based on diffuse hypervoxels as described in any one of claims 1 to 7, characterized in that... include: The image acquisition and preprocessing module is used to receive three-dimensional brain magnetic resonance imaging data, perform preprocessing to obtain processed image data, and use an iterative registration method to register the processed image data to obtain a brain template. The diffusion hypervoxel segmentation module is used to input the brain template into a pre-established water molecule diffusion-based model and output a diffusion-weighted undirected graph. The diffusion geodesic distance is generated based on the diffusion-weighted undirected graph, and the diffusion hypervoxel is segmented based on the diffusion geodesic distance. The fiber bundle clustering module is used to generate a weighted undirected graph of supervoxel neural fiber connections with diffuse supervoxels as vertices, extract the maximum clique of the weighted undirected graph of supervoxel neural fiber connections, and update the weighted undirected graph of supervoxel neural fiber connections based on the maximum clique to achieve fiber bundle clustering. The fiber bundle center curve mapping module is used to register diffuse hypervoxels to individual space, segment individual fiber bundles based on the maximum cluster of individual diffuse hypervoxels, perform curve fitting based on the fiber density image of individual fiber bundles to obtain fiber bundle center curves, extract the features of fiber bundle center curves and perform statistical analysis to realize fiber bundle segmentation analysis.

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 run on a processor. When the processor loads and executes the computer program, it employs the neural fiber tract clustering and segmentation method based on diffuse hypervoxels as described in any one of claims 1 to 7.

10. A computer-readable storage medium storing a computer program, characterized in that, When the computer program is loaded and executed by the processor, it employs the neural fiber bundle clustering and segmentation method based on diffuse hypervoxels as described in any one of claims 1 to 7.

Images