Neural fiber bundle reconstruction method and system based on streamline tracking

Abstract

The invention relates to a nerve fiber bundle reconstruction method and system based on streamline tracking. The method comprises the following steps: extracting a diffusion weighting factor and a diffusion gradient vector of a diffusion weighted image and carrying out preprocessing; constructing a linear equation set for diffusion tensor matrix calculation; solving a diffusion tensor matrix through a non-negatively constrained iterative reweighted least square method; carrying out diagonalization processing on the diffusion tensor matrix, and solving a characteristic value and a characteristic vector; according to the solved characteristic value, calculating the fractional anisotropy of each voxel in the diffusion weighted image, and according to the calculated fractional anisotropy, determining an initial threshold value and an end threshold value of streamline tracking; according to the solved feature vector, the initial threshold value and the termination threshold value, streamline tracking is carried out through an adaptive step length Runge-Kutta integration method; and reconstructing the nerve fiber bundle according to a tracking result. According to the invention, more accurate and efficient nerve fiber bundle reconstruction can be realized.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of medical image processing, and in particular to a neural fiber bundle reconstruction method and system based on streamline tracking. BACKGROUND

[0002] Neural fiber bundle reconstruction based on streamline tracking is based on diffusion tensor imaging (DTI) and other technologies, aiming to visualize the three-dimensional structure of brain white matter fiber bundles, providing key information for understanding brain neural connections. This technology has important value in neurosurgical operation planning, which can provide surgeons with accurate white matter fiber bundle maps to assist in planning safe surgical paths, thereby improving surgical precision and safety and reducing the risk of postoperative neurological damage. In addition, three-dimensional reconstruction of neural fiber bundles is also of great significance for revealing structural changes in white matter fiber bundles, diagnosing and monitoring diseases affecting white matter (such as multiple sclerosis, brain tumor infiltration, etc.), and providing guidance for personalized treatment (such as deep brain stimulation target positioning) of functional neurological diseases (such as Parkinson's disease, epilepsy).

[0003] However, the current mainstream neural fiber bundle reconstruction method based on streamline tracking still faces significant challenges:

[0004] 1) Deterministic tracking algorithm (such as tensor-based streamline tracking): cannot effectively analyze the cross-fiber region, and the tracking error will accumulate along the path, leading to reconstruction trajectory deviation and limited accuracy;

[0005] 2) Probabilistic tracking algorithm (such as methods based on spherical constraint deconvolution): although it can handle cross-fiber problems, the computational complexity is extremely high, and the required time is too long, making it difficult to meet the requirements of efficiency and real-time in clinical practice;

[0006] 3) Deep learning-based algorithm: highly dependent on large-scale and high-quality training data, with high model training cost, and the generalization ability may be restricted by data quality and diversity, and its robustness and universality in clinical application still need to be verified.

[0007] Therefore, it is urgent to develop a neural fiber bundle reconstruction method based on streamline tracking that is computationally efficient, accurate and clinically practical, in order to overcome the shortcomings of existing technologies. SUMMARY

[0008] The present application aims to provide a neural fiber bundle reconstruction method and system based on streamline tracking, which can achieve more accurate and efficient neural fiber bundle reconstruction and has strong clinical practicality.

[0009] In order to achieve the above object, the application provides a neural fiber bundle reconstruction method based on streamline tracking, comprising the following steps:

[0010] S1, extracting a diffusion weighting factor and a diffusion gradient vector of a diffusion weighted image and pre-processing;

[0011] S2, constructing a linear equation set for diffusion tensor matrix calculation;

[0012] S3, solving the diffusion tensor matrix by using an iterative reweighted least square method with non-negative constraint;

[0013] S4, diagonalizing the diffusion tensor matrix to solve eigenvalues and eigenvectors;

[0014] S5, calculating a fractional anisotropy of each voxel in the diffusion weighted image according to the solved eigenvalues, and determining a starting threshold and a terminal threshold for streamline tracking according to the calculated fractional anisotropy;

[0015] S6, performing streamline tracking by using a Runge-Kutta integration method with adaptive step length according to the solved eigenvectors, the starting threshold and the terminal threshold;

[0016] S7, reconstructing a neural fiber bundle according to the tracking result of S6.

[0017] Optionally, in S1, the pre-processing of the extracted diffusion weighting factor and the diffusion gradient vector comprises:

[0018] The diffusion gradient vector is normalized as follows:

[0019]

[0020] wherein g x , g y , and g z represent gradient components of the diffusion gradient vector g in x, y and z directions respectively, ||g max || represents a maximum L2 norm of the gradient direction, and g norm represents a normalized gradient L2 norm;

[0021] According to the principle of diffusion weighted imaging, the square of the diffusion weighting factor and the gradient L2 norm is in linear relationship, and the diffusion weighting factor is normalized as follows:

[0022]

[0023] wherein b o represents an initial diffusion weighting factor.

[0024] Optionally, the S2 comprises:

[0025] According to the Stejskal-Tanner formula, there is

[0026]

[0027] wherein S0 represents a diffusion weighted image without applying any gradient direction, S i represents a diffusion weighted image applying the i-th gradient direction, b i represents a normalized diffusion weighted factor in the i-th gradient direction, D is a diffusion tensor matrix, g i is a gradient vector, representing the i-th gradient direction;

[0028] Taking logarithm on both sides of the above formula and arranging, the calculation formula of the diffusion tensor matrix is as follows:

[0029]

[0030] Optionally, the S3 comprises:

[0031] According to the calculation formula of the diffusion tensor matrix, let the diffusion tensor matrix D be an unknown variable x, be a coefficient matrix A, be y, assuming that there is an unknown variable x, such that Ax-y=0, the solving of the diffusion tensor matrix is converted into an optimization problem as follows:

[0032]

[0033] The objective function of the iterative reweighted least squares method with non-negative constraint is constructed as:

[0034]

[0035] The objective function is expanded as:

[0036]

[0037] wherein, is the i-th row of the coefficient matrix A, y i represents the i-th observation variable, w i is a diagonal element of the weight matrix, and m represents the number of linear equations;

[0038] The initial weight is determined by the signal-to-noise ratio and variance of the diffusion weighted image, that is:

[0039]

[0040] wherein, SNR ia signal-to-noise ratio of an i-th voxel of the diffusion weighted image, representing a variance;

[0041] by the residual The weight w is updated as follows:

[0042]

[0043] where k represents the iteration number, and ∈ represents a minimum value.

[0044] Optionally, the S4 comprises:

[0045] The diffusion tensor matrix is a 3x3 symmetric matrix, so there is an invertible matrix P=[v1, v2, v3] and a diagonal matrix Λ=diag(λ1, λ2, λ3), which satisfy:

[0046] Λ=P -1 DP

[0047] where v1, v2, v3 are three eigenvectors, used to represent the direction of diffusion, and λ1, λ2, λ3 are eigenvalues, used to represent the size of diffusion.

[0048] Optionally, in the S5, the formula for calculating the fractional anisotropy is:

[0049]

[0050] Optionally, the S6 comprises:

[0051] S61, constructing a differential equation based on the obtained eigenvectors, specifically comprising:

[0052] Selecting the eigenvector corresponding to the largest eigenvalue to determine the principal eigenvector of each voxel:

[0053] Based on the principal eigenvector of each voxel point, a vector field is constructed as follows:

[0054] V(r) = V(x, y, z)

[0055] Where r and V(x, y, z) represent the position of each voxel in space and the principal eigenvector;

[0056] Generating streamlines based on the vector field of each voxel point, the streamlines are curves, and the tangent vector of each point on the curve is the same as the direction of the vector field, that is:

[0057]

[0058] Where r(t) represents the position transformation of the streamline, and V(r(t)) represents the direction of the streamline at each point.

[0059] The generation of the streamline is converted into a differential equation solving problem, namely:

[0060]

[0061] Wherein, r0 represents the starting position of the streamline;

[0062] S62, the Runge-Kutta integration method based on adaptive step length is used to solve the differential equation, wherein the definition of the intermediate slope is:

[0063] k1=hr(t n ,r n )

[0064]

[0065] Δ=|y (n,5) -y (n,4) |

[0066] Wherein, t n represents the n th integration step, h represents the integration step length, y (n,4) and y (n,5) represent the approximate solutions of low-order and high-order integrations respectively, and Δ represents the error of the approximate solutions of low-order and high-order integrations;

[0067] S63, the integration step length is dynamically adjusted according to the local error of the current step, which specifically includes:

[0068] The maximum allowed error of the current integration step is selected as:

[0069] e n =|r n |·(1e-6)

[0070] If e n ≤Δ, the current step is accepted, then the integration value of the next step is:

[0071] r n+1 =y (n,5)

[0072] The step length of the next step is:

[0073]

[0074] If e n >Δ, h new is selected as the current step length, and is recalculated until the error meets the condition;

[0075] S64, setting the fractional anisotropy less than or equal to a set threshold as the integration termination condition.

[0076] Optionally, after the streamline reconstruction in the S7, the streamline tracking based neural fiber bundle reconstruction method further comprises:

[0077] Coloring the reconstructed streamlines.

[0078] Based on the same inventive concept, the application further provides a streamline tracking based neural fiber bundle reconstruction system, comprising:

[0079] An extraction module configured to extract diffusion weighting factors and diffusion gradient vectors of diffusion weighted images and perform preprocessing;

[0080] A construction module configured to construct a linear equation set for diffusion tensor matrix calculation;

[0081] A diffusion tensor solving module configured to solve the diffusion tensor matrix by using an iterative reweighted least square method with non-negative constraint;

[0082] A feature solving module configured to perform diagonalization processing on the diffusion tensor matrix to solve eigenvalues and eigenvectors;

[0083] A calculation module configured to calculate fractional anisotropy of each voxel in the diffusion weighted images according to the solved eigenvalues, and determine a starting threshold and a termination threshold for streamline tracking according to the calculated fractional anisotropy;

[0084] A tracking module configured to perform streamline tracking by using a Runge-Kutta integration method with adaptive step length according to the solved eigenvectors, the starting threshold and the termination threshold;

[0085] A reconstruction module configured to perform neural fiber bundle reconstruction and coloring according to the tracking result of the tracking module.

[0086] Based on the same inventive concept, the application further provides a readable storage medium having a computer program stored thereon, wherein the computer program is executable to implement the streamline tracking based neural fiber bundle reconstruction method as described above.

[0087] In the streamline tracking based neural fiber bundle reconstruction method and system provided by the application, at least one of the following beneficial effects is achieved:

[0088] 1) Optimized diffusion tensor calculation: the diffusion tensor is calculated by using an iterative reweighted least square method with non-negative constraint, and the tensor non-negativity is forced and more reasonable weights are given to noise distribution uneven areas, so that the influence of abnormal values is effectively suppressed, and the robustness and accuracy of tensor estimation are improved;

[0089] 2) Adaptive step size integral tracking strategy: the Runge-Kutta integral algorithm with adaptive step size is applied to streamline tracking, and the integral step size is dynamically adjusted through real-time error estimation, which significantly reduces unnecessary calculation under the premise of ensuring tracking accuracy, especially suitable for complex regions with rapid changes in fiber direction (such as cross-fibers), and can effectively maintain the accuracy of the reconstructed path;

[0090] 3) high computational efficiency, reliable accuracy and strong clinical practicability. BRIEF DESCRIPTION OF DRAWINGS

[0091] Those skilled in the art will understand that the provided drawings are used to better understand the present application, and do not constitute any limitation on the scope of the present application. Among them:

[0092] Figure 1 The flowchart of the neural fiber bundle reconstruction method based on streamline tracking provided by an embodiment of the present application;

[0093] Figure 2 The fractional anisotropy map according to the feature vector coloring provided by an embodiment of the present application;

[0094] Figure 3 The left view of the reconstructed neural fiber bundle provided by an embodiment of the present application;

[0095] Figure 4 The top view of the reconstructed neural fiber bundle provided by an embodiment of the present application. DETAILED DESCRIPTION

[0096] In order to make the purpose, advantages and characteristics of the present application more clear, the following will further detail the present application combined with the drawings and specific embodiments. It should be noted that the drawings are greatly simplified and all use non-precise proportions, only to facilitate and clearly assist the purpose of explaining the embodiments of the present application. In order to make the purpose, characteristics and advantages of the present application more obvious and easy to understand, please refer to the drawings. It should be noted that the structure, proportion, size, etc. shown in the drawings attached to the present specification are only used to cooperate with the content disclosed in the specification, so that those skilled in the art can understand and read, and are not used to limit the conditions of the present application. Any modification of structure, change of proportion relationship or adjustment of size, as long as it is the same or similar to the effect and purpose that can be achieved by the present application, should still fall within the scope of the technical content disclosed by the present application.

[0097] Reference Figure 1 The present embodiment provides a neural fiber bundle reconstruction method based on streamline tracking, comprising:

[0098] S1, extracting the diffusion weighting factor and diffusion gradient vector of the diffusion weighted image and pre-processing;

[0099] S2, construct a linear equation set for diffusion tensor matrix calculation;

[0100] S3, solve the diffusion tensor matrix by the non-negative constrained iterative reweighted least squares method;

[0101] S4, diagonalize the diffusion tensor matrix to solve the eigenvalues and eigenvectors;

[0102] S5, according to the solved eigenvalues, calculate the fractional anisotropy of each voxel in the diffusion weighted image, and determine the starting threshold and the ending threshold of the streamline tracking according to the calculated fractional anisotropy;

[0103] S6, according to the solved eigenvectors, the starting threshold and the ending threshold, perform streamline tracking by the Runge-Kutta integration method with adaptive step length;

[0104] S7, according to the tracking result of S6, reconstruct the nerve fiber bundle.

[0105] Specifically, first, S1 is performed, the diffusion weighting factor and the diffusion gradient vector of the diffusion weighted image are extracted and preprocessed. In the present embodiment, when the diffusion tensor estimation is performed based on the diffusion weighted imaging (Diffusion Weighted Imaging, DWI) data, the accurate extraction and standardization of the diffusion gradient parameters (b value and gradient direction) is a key preprocessing step. However, the diffusion gradient parameters of the DWI data generated by different manufacturers' devices are usually stored in different private tags (Private Tags) or standard tags of the DICOM format file, and there is significant format heterogeneity. This inconsistency directly increases the complexity of the diffusion tensor estimation process and may cause processing errors.

[0106] To solve the above problems, the present embodiment provides a unified interface scheme based on the abstract factory design pattern: on the one hand, a standardized method interface is designed, which is used to obtain the normalized diffusion weighting factor (b value) and the diffusion gradient vector (diffusion weighting direction); on the other hand, specific adapter classes (Concrete Adapter) are developed for the DICOM tag storage rules specific to different manufacturers (such as Siemens, GE, Philips, UIH, etc.). These adapters automatically identify the image source manufacturer and accurately analyze and extract the original diffusion weighting factor and diffusion gradient vector data according to the inherent tag specifications of the manufacturer.

[0107] In addition, during the DWI data acquisition process, the gradient vector not only represents the direction, but also may have different intensities (amplitudes). Therefore, by normalizing the gradient vector, the size of the gradient vector in all directions is standardized to 1, thereby eliminating the influence of amplitude difference.

[0108] In this embodiment, the preprocessing of the extracted diffusion weighting factor and the diffusion gradient vector includes:

[0109] The diffuse gradient vector is normalized as follows:

[0110]

[0111] Among them, g x ,g y g z These represent the gradient components of the diffuse gradient vector g in the x, y, and z directions, respectively. max || represents the maximum L2 norm of the gradient direction, g norm This represents the L2 norm of the normalized gradient;

[0112] Based on the principle of diffusion-weighted imaging, the diffusion weighting factor has a linear relationship with the square of the gradient L2 norm. The diffusion weighting factor is normalized as follows:

[0113]

[0114] Among them, b o This represents the initial diffusion weighting factor.

[0115] Then, S2 is executed to construct a system of linear equations for calculating the diffusion tensor matrix. The diffusion tensor reflects the diffusion characteristics of molecules in biological tissues in various directions. Based on the principle of diffusion-weighted imaging, we assume that the diffusion-weighted image without any gradient direction is S0, and the diffusion-weighted image with a certain gradient direction is S... i The diffusion weighting factor is b. i If the diffuse tensor matrix is ​​D, then according to the Stejskal-Tanner formula, we have:

[0116]

[0117] Among them, g i Let be the gradient vector, representing the direction of the i-th gradient;

[0118] Taking the logarithm of both sides of the above equation and rearranging, we obtain the following formula for calculating the diffusion tensor matrix:

[0119]

[0120] According to the above equation, solving the diffusion tensor matrix D can be transformed into a problem of a system of linear equations. The diffusion tensor matrix D is a 3×3 symmetric matrix. In order to solve this equation, it is only necessary to construct a system of six linear equations with at least six non-collinear gradient vectors.

[0121] Then, S3 is performed to solve the diffusion tensor matrix by the non-negative constrained iterative reweighted least squares method. The non-negative constrained iterative reweighted least squares method is used to calculate the diffusion tensor, which can effectively suppress the influence of outliers and improve the robustness and accuracy of tensor estimation by forcing the tensor to be non-negative and giving more reasonable weights to regions with uneven noise distribution.

[0122] In this embodiment, according to the calculation formula of the diffusion tensor matrix, let the diffusion tensor matrix D be an unknown variable x, is a coefficient matrix A, is y, assuming that there is an unknown variable x, such that Ax-y=0, the solution of the diffusion tensor matrix is converted into an optimization problem as follows:

[0123]

[0124] In order to constrain the tensor to be non-negative and give more reasonable weights to regions with uneven noise distribution, effectively suppress the influence of outliers, the objective function of the non-negative constrained iterative reweighted least squares method is constructed as follows:

[0125]

[0126] The objective function is expanded as follows:

[0127]

[0128] wherein, is the i-th row of the coefficient matrix A, y i represents the i-th observation variable, w i is the diagonal element of the weight matrix, and m represents the number of linear equations;

[0129] The initial weight is determined by the signal-to-noise ratio and variance of the diffusion weighted image, that is:

[0130]

[0131] wherein, SNR i represents the signal-to-noise ratio of the i-th voxel of the diffusion weighted image, represents the variance;

[0132] The residual error is updated as follows: The weight w is updated as follows:

[0133]

[0134] wherein, k represents the number of iterations, ∈ represents the minimum value, and is generally taken as 1e-4.

[0135] Then, diagonalization of the diffusion tensor matrix is performed S4 to solve the eigenvalues and eigenvectors. There are many parameters for characterizing the diffusion tensor matrix, but generally, the diffusion tensor matrix is subjected to diagonalization operation, and the eigenvalues represent the size of the diffusion, and the eigenvectors represent the direction of the diffusion. According to the diffusion tensor matrix solved in S3, the eigenvalues and eigenvectors are solved.

[0136] Since the diffusion tensor matrix is a 3x3 symmetric matrix, there is an invertible matrix P = [v1, v2, v3] and a diagonal matrix Λ = diag(λ1, λ2, λ3) that satisfy:

[0137] Λ = P -1 DP

[0138] where v1, v2, v3 are three eigenvectors for representing the direction of diffusion, and λ1, λ2, λ3 are eigenvalues for representing the size of diffusion. Generally, the maximum eigenvalue is taken as the principal eigenvalue, and the corresponding eigenvector represents the main diffusion direction of the fiber bundle.

[0139] Then, S5 is performed to calculate the fractional anisotropy of each voxel in the diffusion-weighted image according to the solved eigenvalues, and to determine the starting threshold and the ending threshold of streamline tracking according to the calculated fractional anisotropy. The fractional anisotropy represents the diffusion degree of the voxel, and the starting threshold and the ending threshold of streamline tracking can be determined according to the fractional anisotropy, which is equivalent to determining the starting threshold and the ending threshold of the subsequent S6 numerical integration based on the adaptive step Runge-Kutta integration method.

[0140] According to the eigenvalues obtained in S4, the fractional anisotropy calculation formula is:

[0141]

[0142] Then, S6 is performed to perform streamline tracking by the adaptive step Runge-Kutta integration method according to the solved eigenvectors, the starting threshold and the ending threshold.

[0143] In this embodiment, S6 includes:

[0144] S61, constructing a differential equation based on the solved eigenvectors, specifically including:

[0145] According to the eigenvectors solved in S4, the eigenvector of each voxel represents the diffusion direction of the current position of the neural fiber, and the eigenvector corresponding to the maximum eigenvalue is selected to determine the principal eigenvector of each voxel:

[0146] Based on the principal eigenvector of each voxel, a vector field can be constructed as follows:

[0147] V(r) = V(x, y, z)

[0148] wherein r and V(x, y, z) represent the position of each voxel in space and the principal eigenvector;

[0149] A streamline is generated based on the vector field of each voxel point, the streamline is a curve, and the tangent vector of each point on the curve is the same as the direction of the vector field, that is:

[0150]

[0151] wherein r(t) represents the position transformation of the streamline, and V(r(t)) represents the direction of the streamline at each point;

[0152] The generation of the streamline is converted into a differential equation solving problem, that is:

[0153]

[0154] wherein r0represents the starting position of the streamline;

[0155] S62, the Runge-Kutta integration method based on adaptive step length is used to solve the differential equation, the Runge-Kutta integration uses multiple intermediate slopes to replace the single movement in the first-order integration, and the accuracy is higher. The Runge-Kutta integration method based on adaptive step length is an effective way to solve the differential equation, and the integration step length is dynamically adjusted according to the local error of the current step, small step length is used in the area with rapid change to ensure the accuracy, and large step length is used in the area with slow change to reduce the calculation amount.

[0156] The intermediate slope is defined as:

[0157] k1 = hr(t n ,r n )

[0158]

[0159] Δ = |y (n,5) -y (n,4) |

[0160] wherein t n represents the nth integration step, h represents the integration step length, y (n,4) and y (n,5) represent the approximate solutions of low-order and high-order integrations respectively, and Δ represents the error of the approximate solutions of low-order and high-order integrations;

[0161] S63, the integration step length is dynamically adjusted according to the local error of the current step, and specifically includes:

[0162] The maximum allowed error of the current integration step is selected as:

[0163] e n = |r n |·(1e-6)

[0164] If e n ≤ Δ, accept the current step, then the next step integral value is:

[0165] r n+1 = y (n,5)

[0166] The next step size is:

[0167]

[0168] If e n > Δ, select h new as the current step size, recalculate until the error meets the condition;

[0169] S64, set the fractional anisotropy less than or equal to a set threshold as the integral termination condition. Generally, select the fractional anisotropy FA≤0.2 as the integral termination condition, at this time it is considered that the current region may be in the brain gray matter region, the neural fiber bundle is less, or the principal eigenvector of the current region is dominated by noise, lacking directionality.

[0170] The embodiment adopts the Runge-Kutta integral algorithm with adaptive step size to perform streamline tracking, and utilizes real-time error estimation to dynamically adjust the integral step size, which significantly reduces unnecessary calculation amount on the premise of ensuring tracking accuracy, and is particularly suitable for complex regions with rapidly changing fiber directions (such as crossing fibers), and can effectively maintain the accuracy of the reconstructed path.

[0171] Finally, S7 is performed, and streamline tracking is performed by the Runge-Kutta integral method with adaptive step size according to the characteristic vector obtained by solving, the starting threshold and the termination threshold. In the embodiment, a seed point r(0) is given as the initial point of numerical integration, and a specific streamline can be generated according to the numerical integration algorithm described in S6. In fiber bundle reconstruction, multiple seed points are generally selected to generate streamline bundles in a specific region to represent the fiber bundle in the specific region.

[0172] In order to visualize the direction distribution of the fiber bundle, after the streamline reconstruction, the neural fiber bundle reconstruction method further includes: coloring the reconstructed streamline. The region along the x direction of the image coordinate system of each fiber bundle is red, the region along the y direction is green, and the region along the z direction is blue. The fractional anisotropy map based on the characteristic vector coloring is shown in Figure 2 , and the reconstructed and colored neural fiber bundle is shown in Figure 3 (left view) and Figure 4 (top view).

[0173] Based on the same inventive concept, the embodiment of the present application further provides a neural fiber bundle reconstruction system based on streamline tracking, comprising:

[0174] an extraction module configured to extract a diffusion weighting factor and a diffusion gradient vector of a diffusion weighted image and perform preprocessing;

[0175] a construction module configured to construct a linear equation set for diffusion tensor matrix calculation;

[0176] a diffusion tensor solving module configured to solve the diffusion tensor matrix by using an iterative reweighted least square method with non-negative constraint;

[0177] a feature solving module configured to perform diagonalization processing on the diffusion tensor matrix to solve eigenvalues and eigenvectors;

[0178] a calculation module configured to calculate fractional anisotropy of each voxel in the diffusion weighted image according to the solved eigenvalues, and determine a starting threshold and a termination threshold for streamline tracking according to the calculated fractional anisotropy;

[0179] a tracking module configured to perform streamline tracking by using a Runge-Kutta integration method with adaptive step length according to the solved eigenvectors, the starting threshold and the termination threshold;

[0180] a reconstruction module configured to perform reconstruction and coloring of neural fiber bundles according to a tracking result of the tracking module.

[0181] Since the neural fiber bundle reconstruction system based on streamline tracking provided by the embodiment of the present application and the neural fiber bundle reconstruction method based on streamline tracking provided above belong to the same inventive concept, the neural fiber bundle reconstruction system based on streamline tracking provided by the present application has all the advantages of the neural fiber bundle reconstruction method based on streamline tracking provided above, and thus the beneficial effects of the neural fiber bundle reconstruction system based on streamline tracking provided by the present application will not be repeated here.

[0182] Based on the same inventive concept, the embodiment of the present application further provides a readable storage medium having a computer program stored thereon, and the computer program can implement the neural fiber bundle reconstruction method based on streamline tracking when executed.

[0183] A readable storage medium can be a tangible device capable of holding and storing instructions for use by an instruction execution device, such as, but not limited to, electrical storage devices, magnetic storage devices, optical storage devices, electromagnetic storage devices, semiconductor storage devices, or any suitable combination thereof. More specific examples of readable storage media (a non-exhaustive list) include: portable computer disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), static random access memory (SRAM), portable compact disc read-only memory (CD-ROM), digital multifunction disc (DVD), memory sticks, floppy disks, mechanical encoding devices, such as punch cards or recessed protrusions storing instructions thereon, and any suitable combination thereof. The computer programs described herein can be downloaded from the readable storage medium to various computing / processing devices, or downloaded via a network, such as the Internet, local area network, wide area network, and / or wireless network, to an external computer or external storage device. Networks can include copper transmission cables, fiber optic transmissions, wireless transmissions, routers, firewalls, switches, gateway computers, and / or edge servers. Each computing / processing device's network adapter card or network interface receives and forwards a computer program from the network for storage on a readable storage medium within the various computing / processing devices. The computer program used to perform the operations of this invention can be assembly instructions, instruction set architecture (ISA) instructions, machine instructions, machine-dependent instructions, microcode, firmware instructions, status setting data, or source code or object code written in any combination of one or more programming languages, including object-oriented programming languages ​​such as Smalltalk, C++, etc., and conventional procedural programming languages ​​such as "C" or similar languages. The computer program can execute entirely on the user's computer, partially on the user's computer, as a standalone software package, partially on the user's computer and partially on a remote computer, or entirely on a remote computer or server. In cases involving remote computers, the remote computer can be connected to the user's computer via any type of network, including a local area network (LAN) or a wide area network (WAN), or it can be connected to an external computer (e.g., via the Internet using an Internet service provider). In some embodiments, electronic circuits, such as programmable logic circuits, field-programmable gate arrays (FPGAs), or programmable logic arrays (PLAs), are personalized by utilizing state information from a computer program. These electronic circuits can execute computer-readable program instructions, thereby realizing various aspects of the invention.

[0184] The computer program can be stored in a tangible computer readable storage medium (or media) having stored thereon instructions that implement the procedures, modules and / or blocks described above. As such, the computer program can be used to implement a computer implemented process that performs the functions / acts specified in the flowchart and / or block diagram.

[0185] The computer program can also be loaded onto a computer, other programmable data processing apparatus, or other device to cause a series of operational steps to be performed on the computer, other programmable apparatus or other device to produce a computer implemented process such that the computer program implemented on the computer, other programmable apparatus, or other device implements the functions / acts specified in the flowchart and / or block diagram.

[0186] Since the readable storage medium provided by the embodiments of the present application belongs to the same inventive concept as the above-mentioned method for reconstructing neural fiber tracts based on streamline tracking, the readable storage medium provided by the present application has all the advantages of the above-mentioned method for reconstructing neural fiber tracts based on streamline tracking, and thus the beneficial effects of the readable storage medium provided by the present application will not be described one by one here.

[0187] In summary, the embodiments of the present application provide a method and system for reconstructing neural fiber tracts based on streamline tracking. On the one hand, the non-negative constrained iterative reweighted least squares method is used to calculate the diffusion tensor, and the tensor non-negativity is forced and a more reasonable weight is given to the uneven noise distribution area, so as to effectively suppress the influence of abnormal values and improve the robustness and accuracy of tensor estimation. On the other hand, the Runge-Kutta integration algorithm with adaptive step size is applied to streamline tracking, the integration step size is dynamically adjusted through real-time error estimation, unnecessary calculation amount is significantly reduced on the premise of ensuring tracking accuracy, and the algorithm is particularly suitable for complex regions with rapid changes in fiber direction (such as crossing fibers), and can effectively maintain the accuracy of the reconstructed path.

[0188] Furthermore, it should be appreciated that the above-described embodiments are merely exemplary as to the present application and need not be practiced as described. It should also be understood that where the above-described embodiments are, for clarity, described in terms of steps or means, these steps or means need not be performed in the order given or even one at a time. Rather, certain steps or means can be performed in an order different than other steps or means or one following another, or performed at least at partially in parallel, to express the sequential order of individual steps or features or functions. Also, description as to the present application need not be performed in the order given or even one at a time. Rather, certain steps or means can be performed in an order different than other steps or means or one following another, or performed at least at partially in parallel, to express the sequential order of individual steps or features or functions. Moreover, it should be recognized that the above-described embodiments need not be performed in the order given or even one at a time. Rather, certain steps or means can be performed in an order different than other steps or means or one following another, or performed at least at partially in parallel, to express the sequential order of individual steps or features or functions.

Claims

1. A method of stream line tracking based neural fiber tract reconstruction, characterized by, The method comprises the following steps: S1, extracting a diffusion weighting factor and a diffusion gradient vector of a diffusion weighted image and preprocessing; S2, constructing a linear equation set for diffusion tensor matrix calculation; S3, solving the diffusion tensor matrix by using an iterative reweighted least squares method with non-negative constraint; S4, diagonalizing the diffusion tensor matrix to solve eigenvalues and eigenvectors; S5, calculating fractional anisotropy of each voxel in the diffusion weighted image according to the solved eigenvalues, and determining a starting threshold and a termination threshold for streamline tracking according to the calculated fractional anisotropy; S6, performing streamline tracking by using a Runge-Kutta integration method with adaptive step length according to the solved eigenvectors, the starting threshold and the termination threshold; S7, reconstructing a nerve fiber bundle according to the tracking result of S6.

2. The stream line tracking based neural fiber bundle reconstruction method of claim 1, wherein, In S1, the preprocessing of the extracted diffusion weighting factor and the diffusion gradient vector comprises: The diffusion gradient vector is normalized as follows: where g x , g y , g z represent the gradient components of the diffusion gradient vector g in the x, y, z directions, respectively, ||g max || represents the maximum L2 norm of the gradient direction, g norm represents the normalized gradient L2 norm; According to the principle of diffusion weighted imaging, the diffusion weighting factor is linearly related to the square of the gradient L2 norm, and the diffusion weighting factor is normalized as follows: where b o represents an initial diffusion weighting factor.

3. The stream line tracking based neural fiber bundle reconstruction method of claim 1, wherein, S2 comprises: According to the Stejskal-Tanner formula, there is: where S0represents a diffusion weighted image without applying any gradient direction, S i represents a diffusion weighted image applying the i-th gradient direction, b i represents a normalized diffusion weighted factor in the i-th gradient direction, D is a diffusion tensor matrix, g i is a gradient vector, representing the i-th gradient direction; Taking the logarithm of both sides of the above formula and arranging it, the calculation formula of the diffusion tensor matrix is as follows:

4. The stream line tracking based neural fiber bundle reconstruction method of claim 3, wherein, S3 comprises: According to the calculation formula of the diffusion tensor matrix, assuming the diffusion tensor matrix D as an unknown variable x, is a coefficient matrix A, is y, assuming that there is an unknown variable x, such that Ax-y=0, the solution of the diffusion tensor matrix is converted into an optimization problem as follows: The objective function of the iterative reweighted least squares method with non-negative constraint is: The objective function is expanded as: wherein is the i-th row of the coefficient matrix A, y i denotes the i-th observed variable, w i is the diagonal element of the weight matrix, m denotes the number of linear equations; The initial weight is determined by the signal-to-noise ratio and the variance of the diffusion weighted image, that is: where SNR i denotes the signal-to-noise ratio of the i-th voxel of the diffusion weighted image, denotes the variance; By residual The weights w are updated as follows: Where k represents the number of iterations, and ∈ represents the minimum value.

5. The stream line tracking based neural fiber bundle reconstruction method of claim 1, wherein, S4 comprises: The diffusion tensor matrix is a 3*3 symmetric matrix, so there is an invertible matrix P=[v1,v2,v3] and a diagonal matrix Λ=diag(λ1,λ2,λ3) that satisfy: Λ = P -1 DP Where v1,v2,v3 are three eigenvectors, which represent the direction of diffusion, and λ1,λ2,λ3 are eigenvalues, which represent the size of diffusion.

6. The stream line tracking based neural fiber bundle reconstruction method of claim 5, wherein, In S5, the fractional anisotropy calculation formula is:

7. The stream line tracking based neural fiber bundle reconstruction method of claim 1, wherein, S6 comprises: S61, constructing a differential equation based on the solved eigenvectors, specifically comprising: Selecting the eigenvector corresponding to the maximum eigenvalue to determine the principal eigenvector of each voxel: Based on the principal eigenvector of each voxel point, a vector field is constructed as follows: V(r) = V(x,y,z) Where r and V(x,y,z) represent the position of each voxel in space and the principal eigenvector; Based on the vector field of each voxel point, a streamline is generated, which is a curve, and the tangent vector of each point on the curve is the same as the direction of the vector field, that is: Where r(t) represents the position transformation of the streamline, and V(r(t)) represents the direction of the streamline at each point; The generation of the streamline is converted into a solution problem of a differential equation, that is: Where r0 represents the starting position of the streamline; S62, solving the differential equation based on the Runge-Kutta integration method with adaptive step length, wherein the definition of the intermediate slope is as follows: k1 = hr(t n ,r n ) Δ = |y (n,5) - y (n,4) | where t n represents the nth integration step, h denotes the integration step size, y (n,4) and y (n,5) represent the low-order and high-order integration approximate solutions, respectively, and Δ denotes the error between the low-order and high-order integration approximate solutions. S63, dynamically adjusting the integration step length according to the local error estimation of the current step, specifically comprising: The maximum allowed error of the current integration step is selected as: e n =|r n |·(1e-6) If e n ≤ Δ, accept the current step, then the integral value of the next step is: r n+1 ＝y (n,5) The step size of the next step is: If e n > Δ, then select h new As the current step size, recalculate until the error satisfies the condition; S64, set the fractional anisotropy less than or equal to the set threshold as the integration termination condition.

8. The stream line tracking based neural fiber bundle reconstruction method of claim 1, wherein, In the S7, after streamline reconstruction, the neural fiber bundle reconstruction method based on streamline tracking further comprises: Coloring the reconstructed neural fiber bundle.

9. A streamline tracking based neural fiber tract reconstruction system, characterized by, Comprise: The extraction module is used for extracting the diffusion weighting factor and the diffusion gradient vector of the diffusion weighted image and pre-processing; The construction module is used for constructing the linear equation set for diffusion tensor matrix calculation; The diffusion tensor solving module is used for solving the diffusion tensor matrix by the non-negative constrained iterative reweighted least square method; The feature solving module is used for diagonalizing the diffusion tensor matrix to solve the eigenvalue and the eigenvector; The calculation module is used for calculating the fractional anisotropy of each voxel in the diffusion weighted image according to the solved eigenvalue, and determining the starting threshold and the termination threshold of the streamline tracking according to the calculated fractional anisotropy; The tracking module is used for performing the streamline tracking by the Runge-Kutta integration method with adaptive step size according to the solved eigenvector, the starting threshold and the termination threshold; The reconstruction module is used for reconstructing and coloring the neural fiber bundle according to the tracking result of the tracking module.

10. A readable storage medium, having stored thereon a computer program, characterized in that, The computer program is executed to realize the neural fiber bundle reconstruction method based on streamline tracking according to any one of claims 1-8. The computer program is executed to realize the neural fiber bundle reconstruction method based on streamline tracking according to any one of claims 1-8.

Images