Magnetic susceptibility tensor imaging method and system with separated microstructure-related frequency shift

Abstract

The present invention discloses a magnetic susceptibility tensor imaging method and system that separates microstructure-related frequency shifts, including: using a three-dimensional flow-compensated gradient echo sequence to scan and collect 16-channel multi-echo complex data S in six head directions, and preprocessing the 16-channel multi-echo complex data S collected in each head direction to obtain multi-echo phase data in six head directions and multi-echo amplitude data in six directions; the six head directions are specifically head anteroposterior, head posterior, head left, head right, head prone, and head supine; based on the multi-echo amplitude data in the head anteroposterior direction, a template M of the brain region of interest is obtained, and the multi-echo amplitude data in other directions are aligned with the echo amplitude data in the head anteroposterior direction to obtain a registration matrix. The method and system of the present application can better describe the magnetic susceptibility distribution of white matter regions and improve the accuracy of quantitative measurement of magnetic susceptibility.

Description

Technical Field

[0001] The present invention relates to the technical field of medical image processing, and in particular to a magnetic susceptibility tensor imaging method and system for separating microstructure-related frequency shifts. Background Art

[0002] Magnetic susceptibility is a fundamental property of matter, reflecting its degree of magnetization in an external magnetic field, and is commonly expressed using the symbol . Hemosiderin, myelin, and deoxyhemoglobin in biological tissues can cause changes in tissue magnetic susceptibility, and certain diseases can also alter tissue magnetic susceptibility. Therefore, quantitatively analyzing the distribution of magnetic susceptibility provides important information for studying tissue function and characteristics, as well as for disease diagnosis and treatment.

[0003] The negative magnetic susceptibility of brain white matter originates from the myelin sheath, which has a multi-membrane structure composed of lipids and proteins, with myelin water filling the spaces between the lipid bilayers. The magnetic susceptibility of a molecule is derived from the sum of the three-dimensional tensors contributed by each molecular bond. Because myelin and proteins are asymmetric molecules, the magnetic susceptibility of myelin exhibits anisotropy. Susceptibility tensor imaging (STI) uses a second-order tensor model to characterize the anisotropy of the magnetic susceptibility of tissues, particularly white matter, and can complement diffusion tensor imaging (DTI) to image white matter fiber bundles in the human brain. STI reconstructs white matter fiber pathways and quantifies changes in white matter myelin content at millimeter-level resolution, obtaining high signal-to-noise ratio magnetic susceptibility imaging.

[0004] The volume magnetic susceptibility of biological tissue can lead to changes in the local magnetic field, thereby producing changes in the resonant frequency, i.e., frequency shifts, within the tissue. Using the phase data of gradient echo imaging, especially at high fields, such changes in the magnetic field or frequency can be detected and high spatial resolution imaging can be performed. In addition to the volume magnetic susceptibility effect, the field perturbations and relaxation rates of the microstructure of white matter and the multiple partitions of white matter fibers (i.e., axons, myelin sheaths, and extracellular spaces) can affect the amplitude and phase of the gradient echo, resulting in non-single exponential attenuation of the gradient echo amplitude and nonlinear changes in the phase. As shown in the attached figure, Figure 1 As shown in Figure 1, within a voxel, a hollow fiber column model is used to simulate the white matter fiber population. The microstructure of the white matter can be represented by two intersecting sets of hollow fiber columns. Different intersecting angles of the fiber columns produce significantly different frequency shifts at the same main magnetic field, H0. This microstructural disturbance in the field pattern cannot be explained by the volume susceptibility tensor model of STI. If not eliminated, it will produce artifacts in the STI image and lead to inaccurate quantitative measurements of magnetic susceptibility. Summary of the Invention

[0005] In order to solve the technical problems existing in the background technology, the present invention proposes a magnetic susceptibility tensor imaging method and system that separates the microstructure-related frequency shift.

[0006] The magnetic susceptibility tensor imaging method proposed in the present invention, which separates the microstructure-related frequency shift, comprises the following steps:

[0007] S1. Use a three-dimensional flow-compensated gradient echo sequence to scan and acquire 16-channel multi-echo complex data S in six head directions, and preprocess the 16-channel multi-echo complex data S acquired in each head direction to obtain multi-echo phase data in six head directions and multi-echo amplitude data in six directions; the six head directions are specifically head anteroposterior, head posterior, head lateral, head lateral, and head supine positions;

[0008] S2. Based on the multi-echo amplitude data in the head anteroposterior direction, a template M of the brain region of interest is obtained, and the multi-echo amplitude data in other directions are registered with the echo amplitude data in the head anteroposterior direction as a reference to obtain a registration matrix;

[0009] S3. Process the multi-echo complex data S collected from the six head directions based on the registration matrix to obtain the multi-echo phase data and amplitude data of the six head directions after registration, and perform unwrapping on the multi-echo phase data from the six head directions after registration based on a fast open source minimum spanning tree algorithm combined with the template package ROMEO template to obtain unwrapped phase data of the multiple echoes in the six head directions;

[0010] S4. Calculate the frequency maps of the multiple echoes in the six head directions based on the phase data of the multiple echoes in the six head directions, and perform background field removal on the frequency maps of the multiple echoes in the six head directions based on the Laplace boundary value background field removal algorithm LBV and the template M of the brain region of interest to obtain the local field maps f of the multiple echoes in the six head directions. local_n ;

[0011] S5, multiple echo local field maps for six head directions local_n Perform frequency fitting based on echo time dependence to obtain local field maps C in six head directions f ;

[0012] S6. Local field maps for six head directions C f Magnetic susceptibility tensor image reconstruction is performed to obtain a magnetic susceptibility tensor distribution map of the head.

[0013] Preferably, the preprocessing specifically includes:

[0014] Taking the absolute value of the 16-channel multi-echo complex data S to obtain a 16-channel multi-echo amplitude signal, and taking the angle of the 16-channel multi-echo complex data S to obtain a 16-channel multi-echo phase signal;

[0015] The square sum of the 16-channel multi-echo amplitude signals is taken and then the square root is taken to obtain the multi-echo amplitude data in six directions of multi-channel combination. The 16-channel multi-echo phase signals are combined using the measured three-dimensional phase offset multi-channel phase combination algorithm MCPC-3D-S to obtain the multi-echo phase data in six directions.

[0016] Preferably, step S2 specifically includes:

[0017] Based on the multi-echo amplitude data in the head anteroposterior direction, the brain extraction toolkit FSL BET of the functional magnetic resonance imaging center software library is used to remove the brain shell of the amplitude map of the first echo in the head anteroposterior direction to obtain the template M of the brain region of interest;

[0018] Taking the third echo amplitude data with good signal-to-noise ratio and contrast in the head anteroposterior direction as the benchmark, the affine transformation toolkit ANTs affine, an automatic medical image registration tool, was used to select the third echo amplitude data in each direction for image registration in other directions to obtain the registration matrix.

[0019] Preferably, step S3 specifically includes:

[0020] Applying the obtained registration matrix to the real and imaginary parts of the six-directional multi-echo complex data S to obtain the registered multi-echo complex data S;

[0021] Taking the absolute value of the registered multi-echo complex data S to obtain the registered multi-echo amplitude data;

[0022] Taking an angle of the registered multi-echo complex data S to obtain the registered multi-echo phase data;

[0023] The registered multi-echo phase data were unwrapped using a fast open-source minimum spanning tree algorithm combined with the ROMEO template method to obtain unwrapped phase data of multiple echoes in six head directions registered to the anteroposterior direction.

[0024] Preferably, the calculating of frequency graphs of the multiple echoes in the six head directions based on the unwrapped phase data of the multiple echoes in the six head directions specifically includes:

[0025] Divide the unwrapped phase data of each echo by 2πTE n To calculate the frequency map of each echo, where TE n is the echo time of n echoes.

[0026] Preferably, step S5 specifically includes:

[0027] Multiple echo local field maps based on six head directions local_nand the echo time TE of the scanning setting, the microstructure-related frequency shift Δf and the local field frequency shift C can be solved using the least squares method. f ;

[0028] Local field frequency shift C f The acquisition process is as follows:

[0029] The effective frequency f of the multi-echo complex data S at the echo time TE is expressed as:

[0030]

[0031] According to the hollow cylindrical fiber model, brain white matter fibers can be divided into three regions: axons, myelin sheaths, and extra-axonal spaces. Therefore, the multi-echo complex data S corresponding to white matter fibers is the sum of the signals of the three regions. The real and imaginary parts of the multi-echo complex data S are expressed as follows:

[0032]

[0033] Among them, v ax ,v m and v ex represents the volume fraction of axon, myelin sheath and extra-axonal space in a voxel, ρ m is the water proton density of the myelin sheath, Δω ax ,Δω m and Δω ex represents the frequency shift caused by the magnetic susceptibility of the axon, myelin sheath, and extra-axonal space, and Represents axons, myelin sheaths, and extraaxonal space time; The simplification is due to the fact that the phase decay effect associated with the magnetic field variation outside the axon is negligible compared to the uniform magnetic field of the axon. Non-myelinated area time;

[0034] Perform small phase approximation and Taylor expansion on sin(x), cos(x) and arctan(x), and only retain the first term of the Taylor series. Formula (1) can be expressed as

[0035]

[0036] Under the condition of short echo time TE, the signal frequency formula (4) can be simplified to

[0037]

[0038] Under long TE conditions, the signal frequency formula (4) can be simplified to

[0039]

[0040] The microstructure-dependent frequency shift Δf can be defined as the difference between the effective frequencies measured in the two states mentioned above:

[0041]

[0042] Combining formulas (4) and (7), the local field pattern signal f of the nth echo is local_n It can be expressed as:

[0043]

[0044] in, and Preferably, step S6 specifically includes:

[0045] Local field map C f The relationship between it and the magnetic susceptibility tensor can be expressed as:

[0046]

[0047] Among them, in m head directions, δB (m) (k) is the k-space form of the local field map obtained by fitting, i, j = 1, 2, 3; m = 1, ... N; Represents the vector of the main magnetic field in the imaging direction coordinate system, k = [k1 k2 k3] is the spatial frequency vector;

[0048] The mathematical expression of the magnetic susceptibility tensor imaging model based on echo time-dependent frequency fitting is:

[0049]

[0050] Among them, the first term is the data fidelity term, and the second term is used to constrain the magnetic susceptibility values ​​of voxels outside the brain area to zero; X ​​is the matrix representation of the magnetic susceptibility tensor:

[0051] F and F- 1 represent Fourier transform and inverse transform operators respectively; M is the template of the brain region of interest;

[0052] Matrix A is the system matrix, expressed as

[0053] FC f Represents the local field map of k-space, FC f =[FC f (1) (k),…FCf (i) (k),…FC f (N) (k)] T , i = 1…N represents the local field map in N directions; λ is the regularization parameter used to adjust the weight of the constraint term;

[0054] The least squares method lsqr is used to solve formula (10), and the iterative convergence condition is set to 1×10 -4 ,The regularization parameter is optimized and set to λ=10 to obtain the magnetic susceptibility tensor distribution map of the head.

[0055] The magnetic susceptibility tensor imaging system proposed in the present invention, which separates the microstructure-related frequency shift, comprises:

[0056] a data scanning and acquisition module for scanning and acquiring 16-channel multi-echo complex data S in six head directions using a three-dimensional flow-compensated gradient echo sequence, and preprocessing the 16-channel multi-echo complex data S acquired in each head direction to obtain multi-echo phase data in the six head directions and multi-echo amplitude data in the six directions; the six head directions are specifically head anteroposterior, head posterior, head lateral, head lateral, and head supine;

[0057] The phase registration module is used to obtain the template M of the brain region of interest based on the multi-echo amplitude data in the head anteroposterior direction, and to register the multi-echo amplitude data in other directions with the echo amplitude data in the head anteroposterior direction as the reference to obtain the registration matrix;

[0058] A dewrapping module is used to process the multi-echo complex data S collected from the six head directions based on the registration matrix to obtain the multi-echo phase data and amplitude data of the six head directions after registration, and to dewrapping the multi-echo phase data of the six head directions after registration based on a fast open source minimum spanning tree algorithm combined with a template package ROMEO template to obtain dewrapped phase data of multiple echoes in the six head directions;

[0059] The background field removal module is used to calculate the frequency map of the multiple echoes in the six head directions based on the phase data of the multiple echoes in the six head directions, and remove the background field of the frequency map of the multiple echoes in the six head directions based on the Laplace boundary value background field removal algorithm LBV and the template M of the brain region of interest to obtain the local field map f of the multiple echoes in the six head directions. local_n ;

[0060] Frequency fitting module for multiple echo local field maps f in six head directions local_n Perform frequency fitting based on echo time dependence to obtain local field maps C in six head directions f ;

[0061] Magnetic susceptibility tensor reconstruction module for local field maps C in six head directions f Magnetic susceptibility tensor image reconstruction is performed to obtain a magnetic susceptibility tensor distribution map of the head.

[0062] The proposed magnetic susceptibility tensor imaging method and system, which separates microstructure-related frequency shifts, separates echo-time-dependent frequency shifts due to microstructure from time-independent frequency shifts due to bulk susceptibility. These shifts are then used for STI reconstruction, enabling more accurate reconstruction and quantification of the magnetic susceptibility tensor. This method can better describe the magnetic susceptibility distribution in white matter regions, including the magnetic susceptibility tensors of single and crossing fiber bundles, thereby improving the accuracy of quantitative susceptibility measurements. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 Schematic diagram of the frequency shift structure generated by the hollow cylindrical fiber model of two groups of fiber groups and fiber groups with different crossing angles in the same main magnetic field;

[0064] Figure 2 Schematic diagram of the STI reconstruction architecture of the magnetic susceptibility tensor imaging method proposed in the present invention that separates the microstructure-related frequency shift;

[0065] Figure 3 This is a schematic diagram of the workflow structure of the magnetic susceptibility tensor imaging method proposed in the present invention that separates the microstructure-related frequency shift;

[0066] Figure 4 Schematic diagram of the system architecture of the magnetic susceptibility tensor imaging system proposed in the present invention that separates the microstructure-related frequency shift. DETAILED DESCRIPTION

[0067] Reference Figure 1-4 The magnetic susceptibility tensor imaging method proposed in the present invention, which separates the microstructure-related frequency shift, comprises the following steps:

[0068] S1. Use a three-dimensional flow-compensated gradient echo sequence to scan and collect 16-channel multi-echo complex data S in six head directions, and preprocess the 16-channel multi-echo complex data S collected in each head direction to obtain six head direction multi-echo phase data and six direction multi-echo amplitude data; the six head directions are specifically head frontal position, head posterior position, head left position, head right position, head prone position, and head supine position.

[0069] In this embodiment, the preprocessing specifically includes:

[0070] Taking the absolute value of the 16-channel multi-echo complex data S to obtain a 16-channel multi-echo amplitude signal, and taking the angle of the 16-channel multi-echo complex data S to obtain a 16-channel multi-echo phase signal;

[0071] The square sum of the 16-channel multi-echo amplitude signals is taken and then the square root is taken to obtain the multi-echo amplitude data in six directions of multi-channel combination. The 16-channel multi-echo phase signals are combined using the measured three-dimensional phase offset multi-channel phase combination algorithm MCPC-3D-S to obtain the multi-echo phase data in six directions.

[0072] S2. Based on the multi-echo amplitude data in the head anteroposterior direction, a template M of the brain region of interest is obtained, and the multi-echo amplitude data in other directions are registered with the echo amplitude data in the head anteroposterior direction as a reference to obtain a registration matrix.

[0073] In this embodiment, step S2 specifically includes:

[0074] Based on the multi-echo amplitude data in the head anteroposterior direction, the brain extraction toolkit FSL BET of the functional magnetic resonance imaging center software library is used to remove the brain shell of the amplitude map of the first echo in the head anteroposterior direction to obtain the template M of the brain region of interest;

[0075] Taking the third echo amplitude data with good signal-to-noise ratio and contrast in the head anteroposterior direction as the benchmark, the affine transformation toolkit ANTs affine, an automatic medical image registration tool, was used to select the third echo amplitude data in each direction for image registration in other directions to obtain the registration matrix.

[0076] S3. Process the multi-echo complex data S collected from the six head directions based on the registration matrix to obtain the multi-echo phase data and amplitude data of the six head directions after registration, and perform dewrapping on the multi-echo phase data of the six head directions after registration based on the fast open source minimum spanning tree algorithm combined with the template package ROMEO template to obtain the dewrapped phase data of multiple echoes in the six head directions.

[0077] In this embodiment, step S3 specifically includes:

[0078] Applying the obtained registration matrix to the real and imaginary parts of the six-directional multi-echo complex data S to obtain the registered multi-echo complex data S;

[0079] Taking the absolute value of the registered multi-echo complex data S to obtain the registered multi-echo amplitude data;

[0080] Taking an angle of the registered multi-echo complex data S to obtain the registered multi-echo phase data;

[0081] The registered multi-echo phase data were unwrapped using a fast open-source minimum spanning tree algorithm combined with the ROMEO template method to obtain unwrapped phase data of multiple echoes in six head directions registered to the anteroposterior direction.

[0082] S4. Calculate the frequency maps of the multiple echoes in the six head directions based on the phase data of the multiple echoes in the six head directions, and perform background field removal on the frequency maps of the multiple echoes in the six head directions based on the Laplace boundary value background field removal algorithm LBV and the template M of the brain region of interest to obtain the local field maps f of the multiple echoes in the six head directions. local_n .

[0083] In this embodiment, frequency maps of the multiple echoes in the six head directions are calculated based on the unwrapped phase data of the multiple echoes in the six head directions, specifically including:

[0084] Divide the unwrapped phase data of each echo by 2πTE n To calculate the frequency map of each echo, where TE n is the echo time of n echoes.

[0085] In this embodiment, the collected data in other directions are aligned to the head position of the subject using the ANTs affine linear alignment toolkit, and the third echo amplitude data with good signal-to-noise ratio and contrast in each direction is selected for alignment. The obtained alignment matrix is ​​applied to the real and imaginary parts of the multi-echo complex data, and then the absolute value and angle of the registered complex data are taken to obtain the amplitude and phase information. The aligned multi-echo phase is unwrapped using the ROMEO template method. The frequency map containing only the local field disturbance is fitted using a frequency fitting method based on echo time dependence to obtain a local field map that removes the microstructure frequency shift and only contains the frequency shift caused by tissue susceptibility, which is input into the susceptibility tensor model.

[0086] S5, multiple echo local field maps for six head directions local_n Perform frequency fitting based on echo time dependence to obtain local field maps C in six head directions f .

[0087] In this embodiment, step S5 specifically includes:

[0088] Multiple echo local field maps based on six head directions local_n and the echo time TE of the scanning setting, the microstructure-related frequency shift Δf and the local field frequency shift C can be solved using the least squares method. f ;

[0089] Local field frequency shift C f The acquisition process is as follows:

[0090] The effective frequency f of the multi-echo complex data S at the echo time TE is expressed as:

[0091]

[0092] According to the hollow cylindrical fiber model, brain white matter fibers can be divided into three regions: axons, myelin sheaths, and extra-axonal spaces. Therefore, the multi-echo complex data S corresponding to white matter fibers is the sum of the signals of the three regions. The real and imaginary parts of the multi-echo complex data S are expressed as follows:

[0093]

[0094] Among them, v ax ,v m and v ex represents the volume fraction of axon, myelin sheath and extra-axonal space in a voxel, ρ m is the water proton density of the myelin sheath, Δω ax ,Δω m and Δω ex represents the frequency shift caused by the magnetic susceptibility of the axon, myelin sheath, and extra-axonal space, and Represents axons, myelin sheaths, and extraaxonal space time; The simplification is due to the fact that the phase decay effect associated with the magnetic field variation outside the axon is negligible compared to the uniform magnetic field of the axon. Non-myelinated area time;

[0095] Perform small phase approximation and Taylor expansion on sin(x), cos(x) and arctan(x), and only retain the first term of the Taylor series. Formula (1) can be expressed as

[0096]

[0097] Under the condition of short echo time TE, the signal frequency formula (4) can be simplified to

[0098]

[0099] Under long TE conditions, the signal frequency formula (4) can be simplified to

[0100]

[0101] The microstructure-dependent frequency shift Δf can be defined as the difference between the effective frequencies measured in the two states mentioned above:

[0102]

[0103] Combining formulas (4) and (7), the local field pattern signal f of the nth echo is local_n It can be expressed as:

[0104]

[0105] in, and

[0106] S6. Local field maps for six head directions C f Magnetic susceptibility tensor image reconstruction is performed to obtain a magnetic susceptibility tensor distribution map of the head.

[0107] In this embodiment, step S6 specifically includes:

[0108] Local field map C f The relationship between it and the magnetic susceptibility tensor can be expressed as:

[0109]

[0110] Among them, in m head directions, δB (m) (k) is the k-space form of the local field map obtained by fitting, i,j=1,2,3; m=1,...N; Represents the vector of the main magnetic field in the imaging direction coordinate system, k = [k1 k2 k3] is the spatial frequency vector;

[0111] The mathematical expression of the magnetic susceptibility tensor imaging model based on echo time-dependent frequency fitting is:

[0112]

[0113] Among them, the first term is the data fidelity term, and the second term is used to constrain the magnetic susceptibility values ​​of voxels outside the brain area to zero; X ​​is the matrix representation of the magnetic susceptibility tensor:

[0114] F and F- 1 represent Fourier transform and inverse transform operators respectively; M is the template of the brain region of interest;

[0115] Matrix A is the system matrix, expressed as

[0116] FC f Represents the local field map of k-space, FC f =[FC f (1) (k),…FCf (i) (k), …FC f (N) (k)] T , i = 1…N represents the local field map in N directions; λ is the regularization parameter used to adjust the weight of the constraint term;

[0117] The least squares method lsqr is used to solve formula (10), and the iterative convergence condition is set to 1×10- 4 ,The regularization parameter is optimized and set to λ=10 to obtain the magnetic susceptibility tensor distribution map of the head.

[0118] In this embodiment, susceptibility tensor imaging (STI) is an emerging magnetic resonance imaging technique used to quantitatively measure and image the anisotropic magnetic susceptibility distribution of tissues. Magnetic susceptibility reflects the magnetization intensity of a substance under the action of an external magnetic field. STI not only provides information about the magnitude of the magnetic susceptibility, but also its direction, and is therefore called "tensor" imaging. For certain anisotropic materials, such as brain white matter fibers, the magnetic susceptibility can be approximated as a real symmetric 3×3 tensor, expressed as χ:

[0119]

[0120] For axisymmetric structures, this tensor contains 6 independent elements: 11 , χ 12 , χ 13 , χ 22 , χ 23 , χ 33 .

[0121] The goal of STI is to characterize the magnetic susceptibility tensor at each voxel in the brain. Brain tissue generates its own magnetic field, which perturbs the main magnetic field. This local magnetic field change causes a change in the resonant frequency, which can be measured using the phase image in the gradient echo. The local magnetic field perturbation measured by the phase image can be obtained by convolving the magnetic susceptibility distribution with the unit dipole. Therefore, the magnetic field perturbation ΔB(r) and the three-dimensional magnetic susceptibility distribution χ can be expressed as:

[0122]

[0123] where ΔB(r) is the frequency shift of the field in Hz, represents the vector of the main magnetic field in the subject's coordinate system, k = [k x k y k z ]'and is the spatial frequency shift vector, F and F -1Represent the Fourier transform and inverse transform operators respectively. So in k-space or frequency domain it can be expressed as:

[0124] δB(k)=a 11 χ 11 (k)+a 12 χ 12 (k)+a 13 χ 13 (k)+a 22 χ 22 (k)+a 23 χ 23 (k)+a 33 χ 33 (k); (11)

[0125] If data from N head directions are collected, the relationship between the local magnetic field disturbance and the magnetic susceptibility tensor can be expanded as follows:

[0126]

[0127] in, i, j=1,2,3; m=1,...N;

[0128] Solving the magnetic susceptibility tensor requires phase data in at least six head directions. Since the main magnetic field of the magnetic resonance imaging device cannot rotate, the subject's brain needs to be rotated at multiple angles inside the head coil to achieve this.

[0129] In this embodiment, the STI reconstruction process specifically includes: first, collecting multi-channel three-dimensional multi-echo gradient echo signals from N (N≥6) head directions, combining the amplitude and phase signals of the multi-channels, and removing the initial phase offset of the echo time TE=0ms from the phase signal. Secondly, aligning the combined amplitude and phase maps of other directions to the subject's head in the frontal position, and then performing phase unwrapping on the phase map, and dividing the obtained phase map by 2πTE n (TE n After removing the background field from the frequency map, the frequency map C in N directions is obtained by frequency fitting based on multi-echo time dependence. f , which is input into the magnetic susceptibility tensor model as a local field map, can solve the magnetic susceptibility tensor and further obtain the mean magnetic susceptibility (MMS) and magnetic susceptibility anisotropy (MSA) of the brain, and the magnetic susceptibility map parallel to the direction of white matter fibers χ || and the magnetic susceptibility map perpendicular to the direction of white matter fibers χ ⊥ .

[0130] In this embodiment, frequency maps can be combined to achieve more accurate field perturbation estimation than single echo frequency maps by combining multiple echo frequencies. This is because combining multiple echo frequency maps can remove the initial phase offset at echo time TE = 0 ms and provide a higher signal-to-noise ratio in reconstructing local field maps and magnetic susceptibility maps. Existing echo combination methods include weighted averaging, linear fitting, and nonlinear fitting.

[0131] Non-linear fitting uses complex data for fitting, taking into account the Gaussian noise in the complex image, and fitting the complex signals of multiple echo times, that is, the exponential expression of the signal, and fitting the frequency shift of the field and the initial phase offset as parameters. This method requires three or more echo data for fitting. Linear fitting is to first solve the phase diagram of each echo complex signal, unwrap it in the time dimension, and then obtain the frequency diagram through weighted least squares fitting. Generally, the amplitude corresponding to different echoes is used as the weight. The weighted averaging method requires spatial expansion of the data of each echo, and subtracting the initial phase offset from the phase of each echo to perform weighted averaging to solve the frequency diagram. The weight is the echo time × amplitude. The above three methods are all based on the theory that the phase changes linearly with the echo time TE, that is The nonlinear evolution of the gradient echo phase caused by the frequency shift caused by the microstructure of the white matter is not considered. If the frequency shift caused by the microstructure is not properly considered in the existing STI model, artifacts will be introduced in the STI reconstruction image, and the tissue susceptibility may be overestimated or underestimated. This modeling error may lead to different susceptibility values ​​observed in different brain regions that are dependent on the echo time, that is, different susceptibility tensor maps reconstructed with different echo times, as well as differences in the reconstruction of susceptibility tensor maps under different protocols, different magnetic resonance scanners and different field strength magnetic fields. Therefore, the present invention considers the microstructure of white matter fibers, and uses formula (11) constructed by a hollow fiber column model of multiple partitions (axons, myelin sheaths and extracellular space) to solve C by the local field frequency f and echo time TE of the multiple echoes through a nonlinear least squares solver. f , removing the frequency shift caused by white matter microstructure Then it is input into the magnetic susceptibility tensor model and the least squares LSQR solution is used to obtain a more accurate magnetic susceptibility tensor reconstructed image.

[0132] Reference Figure 1-4 The magnetic susceptibility tensor imaging system for separating microstructure-related frequency shifts proposed in the present invention includes:

[0133] a data scanning and acquisition module for scanning and acquiring 16-channel multi-echo complex data S in six head directions using a three-dimensional flow-compensated gradient echo sequence, and preprocessing the 16-channel multi-echo complex data S acquired in each head direction to obtain multi-echo phase data in the six head directions and multi-echo amplitude data in the six directions; the six head directions are specifically head anteroposterior, head posterior, head lateral, head lateral, and head supine;

[0134] The phase registration module is used to obtain the template M of the brain region of interest based on the multi-echo amplitude data in the head anteroposterior direction, and to register the multi-echo amplitude data in other directions with the echo amplitude data in the head anteroposterior direction as the reference to obtain the registration matrix;

[0135] A dewrapping module is used to process the multi-echo complex data S collected from the six head directions based on the registration matrix to obtain the multi-echo phase data and amplitude data of the six head directions after registration, and to dewrapping the multi-echo phase data of the six head directions after registration based on a fast open source minimum spanning tree algorithm combined with a template package ROMEO template to obtain dewrapped phase data of multiple echoes in the six head directions;

[0136] The background field removal module is used to calculate the frequency map of the multiple echoes in the six head directions based on the phase data of the multiple echoes in the six head directions, and remove the background field of the frequency map of the multiple echoes in the six head directions based on the Laplace boundary value background field removal algorithm LBV and the template M of the brain region of interest to obtain the local field map f of the multiple echoes in the six head directions. local_n ;

[0137] Frequency fitting module for multiple echo local field maps f in six head directions local_n Perform frequency fitting based on echo time dependence to obtain local field maps C in six head directions f ;

[0138] Magnetic susceptibility tensor reconstruction module for local field maps C in six head directions f Magnetic susceptibility tensor image reconstruction is performed to obtain a magnetic susceptibility tensor distribution map of the head.

[0139] The above description is only a preferred specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any technician familiar with the technical field, within the technical scope disclosed by the present invention, who makes equivalent replacements or changes based on the technical solution and inventive concept of the present invention, should be covered by the scope of protection of the present invention.

Claims

1. A magnetic susceptibility tensor imaging method that separates microstructure-related frequency shifts, characterized in that: The following steps are involved: S1. Use a three-dimensional flow-compensated gradient echo sequence to scan and acquire 16-channel multi-echo complex data S in six head directions, and preprocess the 16-channel multi-echo complex data S acquired in each head direction to obtain multi-echo phase data in six head directions and multi-echo amplitude data in six directions; the six head directions are specifically head anteroposterior, head posterior, head lateral, head lateral, and head supine positions; S2. Based on the multi-echo amplitude data in the head anteroposterior direction, a template M of the brain region of interest is obtained, and the multi-echo amplitude data in other directions are registered with the echo amplitude data in the head anteroposterior direction as a reference to obtain a registration matrix; S3. Process the multi-echo complex data S collected from the six head directions based on the registration matrix to obtain the multi-echo phase data and amplitude data of the six head directions after registration, and perform unwrapping on the multi-echo phase data from the six head directions after registration based on a fast open source minimum spanning tree algorithm combined with the template package ROMEO template to obtain unwrapped phase data of the multiple echoes in the six head directions; S4. Calculate the frequency maps of the multiple echoes in the six head directions based on the phase data of the multiple echoes in the six head directions, and perform background field removal on the frequency maps of the multiple echoes in the six head directions based on the Laplace boundary value background field removal algorithm LBV and the template M of the brain region of interest to obtain the local field maps f of the multiple echoes in the six head directions. local_n ; S5, multiple echo local field maps for six head directions local_n Perform frequency fitting based on echo time dependence to obtain local field maps C in six head directions f ; S6. Local field maps for six head directions C f Performing magnetic susceptibility tensor image reconstruction to obtain a magnetic susceptibility tensor distribution map of the head; Step S5 specifically includes: Multiple echo local field maps based on six head directions local_n and the echo time TE of the scanning setting, and the least squares method is used to solve the microstructure related frequency shift Δf and the local field frequency shift C f ; Local field frequency shift C f The acquisition process is as follows: The effective frequency f of the multi-echo complex data S at the echo time TE is expressed as: According to the hollow cylindrical fiber model, brain white matter fibers can be divided into three regions: axons, myelin sheaths, and extraaxonal spaces. Therefore, the multi-echo complex data S corresponding to white matter fibers is the sum of the signals of the three regions. The real and imaginary parts of the multi-echo complex data S are expressed as follows: Among them, v ax ,v m and v ex represents the volume fraction of axon, myelin sheath and extra-axonal space in a voxel, ρ m is the water proton density of the myelin sheath, Δω ax ,Δω m and Δω ex represents the frequency shift caused by the magnetic susceptibility of the axon, myelin sheath, and extra-axonal space, and Represents axons, myelin sheaths, and extraaxonal space time; The simplification is due to the negligible phase decay effect associated with the magnetic field variation outside the axon compared to the uniform magnetic field of the axon. Non-myelinated area time; Perform small phase approximation and Taylor expansion on sin(x), cos(x) and arctan(x), and only retain the first term of the Taylor series. Formula (1) can be expressed as Under the condition of short echo time TE, the signal frequency formula (4) is simplified to Under long TE conditions, the signal frequency formula (4) is simplified to The microstructure-dependent frequency shift Δf is defined as the difference between the effective frequencies measured in the two states mentioned above: Combining formulas (4) and (7), the local field pattern signal f of the nth echo is local_n Expressed as: in, and TE n Indicates the echo time of the nth echo.

2. The magnetic susceptibility tensor imaging method with separated microstructure-related frequency shifts according to claim 1, characterized in that: The pretreatment specifically includes: Taking the absolute value of the 16-channel multi-echo complex data S to obtain a 16-channel multi-echo amplitude signal, and taking the angle of the 16-channel multi-echo complex data S to obtain a 16-channel multi-echo phase signal; The square sum of the 16-channel multi-echo amplitude signals is taken and then the square root is taken to obtain the multi-echo amplitude data in six directions of multi-channel combination. The 16-channel multi-echo phase signals are combined using the measured three-dimensional phase offset multi-channel phase combination algorithm MCPC-3D-S to obtain the multi-echo phase data in six directions.

3. The magnetic susceptibility tensor imaging method with separated microstructure-related frequency shifts according to claim 1, characterized in that: Step S2 specifically includes: Based on the multi-echo amplitude data in the head anteroposterior direction, the brain extraction toolkit FSL BET of the functional magnetic resonance imaging center software library is used to remove the brain shell of the amplitude map of the first echo in the head anteroposterior direction to obtain the template M of the brain region of interest; Taking the third echo amplitude data with good signal-to-noise ratio and contrast in the head anteroposterior direction as the benchmark, the affine transformation toolkit ANTs affine, an automatic medical image registration tool, was used to select the third echo amplitude data in each direction for image registration in other directions to obtain the registration matrix.

4. The magnetic susceptibility tensor imaging method with separated microstructure-related frequency shifts according to claim 1, characterized in that: Step S3 specifically includes: Applying the obtained registration matrix to the real and imaginary parts of the six-directional multi-echo complex data S to obtain the registered multi-echo complex data S; Taking the absolute value of the registered multi-echo complex data S to obtain the registered multi-echo amplitude data; Taking an angle of the registered multi-echo complex data S to obtain the registered multi-echo phase data; The registered multi-echo phase data were unwrapped using a fast open-source minimum spanning tree algorithm combined with the ROMEO template method to obtain unwrapped phase data of multiple echoes in six head directions registered to the anteroposterior direction.

5. The magnetic susceptibility tensor imaging method with separated microstructure-related frequency shifts according to claim 1, characterized in that: The step of calculating frequency graphs of the multiple echoes in the six head directions based on the unwrapped phase data of the multiple echoes in the six head directions specifically includes: Divide the unwrapped phase data of each echo by 2πTE n To calculate the frequency map of each echo, where TE n is the echo time of n echoes.

6. The magnetic susceptibility tensor imaging method with separated microstructure-related frequency shifts according to claim 1, characterized in that: Step S6 specifically includes: Local field map C f The relationship between and the magnetic susceptibility tensor can be expressed as: Among them, in m head directions, δB (m) (k) is the k-space form of the local field map obtained by fitting, Represents the vector of the main magnetic field in the imaging direction coordinate system, k = [k1 k2 k3] is the spatial frequency vector; The mathematical expression of the magnetic susceptibility tensor imaging model based on echo time-dependent frequency fitting is: Among them, the first term is the data fidelity term, and the second term is used to constrain the magnetic susceptibility values ​​of voxels outside the brain area to zero; X ​​is the matrix representation of the magnetic susceptibility tensor: F and F -1 represent Fourier transform and inverse transform operators respectively; M is the template of the brain region of interest; Matrix A is the system matrix, expressed as FC f Represents the local field map of k-space, FC f =[FC f (1) (k),…FC f (i) (k),…FC f (N) (k)] T ,i=1…N represents the local field map in N directions;λ is the regularization parameter used to adjust the weight of the constraint term; The least squares method lsqr is used to solve formula (10), and the iterative convergence condition is set to 1×10 -4 ,The regularization parameter is optimized and set to λ=10 to obtain the magnetic susceptibility tensor distribution map of the head.

7. A magnetic susceptibility tensor imaging system that separates microstructure-related frequency shifts, characterized in that: The method for magnetic susceptibility tensor imaging with separation of microstructure-related frequency shifts according to any one of claims 1 to 6, wherein the system comprises: a data scanning and acquisition module for scanning and acquiring 16-channel multi-echo complex data S in six head directions using a three-dimensional flow-compensated gradient echo sequence, and preprocessing the 16-channel multi-echo complex data S acquired in each head direction to obtain multi-echo phase data in the six head directions and multi-echo amplitude data in the six directions; the six head directions are specifically head anteroposterior, head posterior, head lateral, head lateral, and head supine; The phase registration module is used to obtain the template M of the brain region of interest based on the multi-echo amplitude data in the head anteroposterior direction, and to register the multi-echo amplitude data in other directions with the echo amplitude data in the head anteroposterior direction as the reference to obtain the registration matrix; A dewrapping module is used to process the multi-echo complex data S collected from the six head directions based on the registration matrix to obtain the multi-echo phase data and amplitude data of the six head directions after registration, and to dewrapping the multi-echo phase data of the six head directions after registration based on a fast open source minimum spanning tree algorithm combined with a template package ROMEO template to obtain dewrapped phase data of multiple echoes in the six head directions; The background field removal module is used to calculate the frequency map of the multiple echoes in the six head directions based on the phase data of the multiple echoes in the six head directions, and remove the background field of the frequency map of the multiple echoes in the six head directions based on the Laplace boundary value background field removal algorithm LBV and the template M of the brain region of interest to obtain the local field map f of the multiple echoes in the six head directions. local_n ; Frequency fitting module for multiple echo local field maps f in six head directions local_n Perform frequency fitting based on echo time dependence to obtain local field maps C in six head directions f ; Magnetic susceptibility tensor reconstruction module for local field maps C in six head directions f Magnetic susceptibility tensor image reconstruction is performed to obtain a magnetic susceptibility tensor distribution map of the head.

Images