Method for calculating quantitative magnetic susceptibility image through total field inversion based on Green function

Abstract

The invention discloses a method for calculating a quantitative magnetic susceptibility image through total field inversion based on a Green function. The method is used for separating a background magnetic field and reconstructing a whole-brain quantitative magnetic susceptibility image in magnetic resonance imaging. According to the boundary element method full-field inversion method, filtering and regularization processing on tissue magnetic susceptibility are prevented from being additionally introduced when background field interference is removed, the background field can be effectively separated, a whole-brain quantitative magnetic susceptibility image can be reconstructed, the boundary erosion phenomenon cannot occur, and the boundary element method full-field inversion method particularly shows excellent performance in a cortex area. According to the method, a fast multipole method is adopted to efficiently process storage and calculation of a dense matrix generated by a boundary integral equation.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of magnetic resonance imaging, and in particular to a method for calculating quantitative susceptibility images based on Green's function total field inversion. BACKGROUND

[0002] Quantitative Susceptibility Mapping (QSM) is an advanced magnetic resonance imaging technique that quantifies the magnetic susceptibility of tissues by deconvolving the signal phase data. In neuroscience research, QSM can identify and quantify sources of magnetic susceptibility in tissues, such as iron, myelin, deoxyhemoglobin iron, and calcification, with wide application prospects. However, in practical applications, background field removal is a crucial step in QSM processing, especially in brain surface areas such as the cerebral cortex, where the interference of the background field is very significant.

[0003] Existing background field removal methods often have errors near the brain tissue boundary, resulting in unsatisfactory removal results. These errors are usually caused by inaccurate prior assumptions about the background field or loss of low-frequency information, affecting the analysis results of the entire brain region and cortex. Background field removal is crucial for minimizing errors in QSM, especially in brain boundary regions, where the magnetic susceptibility difference is significant and susceptible to background field interference. SUMMARY

[0004] The purpose of the present application is to provide a method for calculating quantitative susceptibility images based on Green's function total field inversion, which uses the Green's function of the Laplace equation to establish an integral expression of the background field, and then uses the total field to reconstruct the magnetic susceptibility image. The present application accurately simulates the background field by constructing a surface triangular model and combining the fast multipole expansion technique, avoiding the problems of hypothesis error and loss of low-frequency information in traditional background field removal methods. Compared with existing methods, the present application can effectively improve the image quality of the brain boundary region, reduce the interference of the background field, and make the magnetic susceptibility reconstruction of the cerebral cortex and other key regions more accurate and uniform.

[0005] The technical solution adopted to achieve the above-mentioned purposes of the present application is as follows:

[0006] A method for calculating quantitative susceptibility images based on Green's function total field inversion, characterized in that the method comprises the following steps:

[0007] Step 1: acquire magnetic resonance raw images using a multi-echo gradient echo sequence;

[0008] Step 2: preprocess the modulus and phase images of all echoes to obtain the region of interest and total field images;

[0009] Step 3: Calculate the surface triangulation model based on the region of interest obtained in Step 2;

[0010] Step 4: Using the surface triangular model obtained in Step 3, combined with other sequence-related parameters, construct the optimized formula for the total field inversion, expressed as the following equation:

[0011] (1)

[0012] in, The desired magnetic susceptibility value is... The boundary conditions are to be determined. for intermediate variables, For preprocessing, For magnitude weighting, It is a dipole nucleus. For the integral operator of the Laplace equation, The total field obtained in step 2, It is the regularization weight. It is a structural information mask obtained from the model diagram.

[0013] In formula (1), the operator The following system of integral equations, constructed using Green's function, is obtained:

[0014] (2)

[0015] in, is the Green's function of the Laplace equation, S is the brain surface, and n′ is the unit outward normal vector of S'. In this method, the algorithm uses the brain surface triangulation obtained from step 3 to discretize equation (2) to obtain the operator. Discrete matrix form Used for subsequent numerical calculations. It has the following matrix form:

[0016] (3)

[0017] Formula (3) is The specific composition consists of four matrices: T, Q, T', and Q', with elements being the integrals of the surface triangular model faces. -1 The expression for each element of T, Q, T', Q' is also given in formula (3). The subscripts i and j represent the numbers of the surface triangles, and (x, y, z) are the coordinates in the three-dimensional image space, used to represent the numbers of the matrix elements.

[0018] Step 5: Compress the matrix using the fast multipole expansion algorithm. The size of the distribution map is so small that it can be stored in a desktop workstation or a small server and used for further calculation.

[0019] Step 6: The Gauss-Newton method is used to numerically calculate the solution of formula (1) in step 4 to obtain the magnetic susceptibility distribution map.

[0020] Compared with the prior art, the technical scheme of the present application can achieve the following effects:

[0021] 1) The obtained magnetic susceptibility distribution map is more accurate, and the interference of the background magnetic field is eliminated;

[0022] 2) The magnetic susceptibility values of the brain cortex and white matter and other regions close to the brain boundary obtained by the present application have high accuracy, and the tissue is more uniform;

[0023] 3) The present application is suitable for quantitative measurement of brain iron and other magnetic susceptibility sources in clinical radiology. BRIEF DESCRIPTION OF DRAWINGS

[0024] Figure 1 The flowchart of the present application is shown in the figure;

[0025] Figure 2 The schematic diagram of the magnetic field physical model constructed by the present application is shown in the figure; the inner layer of the schematic diagram represents the magnetic susceptibility , and the white grid wrapped in the outer layer is the triangular surface modeling of the brain surface;

[0026] Figure 3 The magnetic susceptibility distribution map obtained by the present application is shown in the figure. DETAILED DESCRIPTION

[0027] The present application will be further described in detail in combination with the following specific examples and drawings. The process, conditions, experimental scheme and method for implementing the present application are all general knowledge and common sense in the art, and the present application does not have special limitations.

[0028] Reference Figure 1 and Figure 2The total field inversion method for magnetic resonance imaging based on Green function provided by the application improves the accuracy of the magnetization imaging by simulating the background field through the construction of octree structure and the combination of fast multipole expansion technology. Before the total field inversion, the boundary of the region of interest is modeled first, and the background field is calculated accurately by combining the boundary integral method and the Green function model. The inversion process is constrained by introducing the boundary condition, so as to ensure the high-precision reconstruction of the whole brain magnetization without eroding the boundary region of the brain. The method can effectively reduce the error in the traditional background field removal method, and improve the signal-to-noise ratio and imaging quality of the magnetization distribution map of the cerebral cortex and other regions close to the boundary of the brain.

[0029] The following describes the specific implementation process of the background field modeling and total field inversion calculation of the original data obtained by the multi-echo gradient echo sequence acquisition in the application. The magnetic resonance image data is acquired by the multi-echo gradient echo sequence on a 3T magnetic resonance imaging device system (Siemens MAGNETOM Prisma Fit 3T).

[0030] Step 1: The multi-echo gradient echo sequence is used to scan the subject's brain on the 3T magnetic resonance imaging device system to obtain the multi-echo original image, including the phase diagram and the modulus diagram.

[0031] The specific implementation of scanning the subject's brain by the multi-echo gradient echo sequence is the general process of magnetic resonance scanning, and the obtained original image is a plurality of three-dimensional images acquired at different echo times TE. In this embodiment, the magnetic resonance scanning parameters are TR / TE1 / ΔTE = 38 ms / 6.4 ms / 5.2 ms, the number of echoes is 6, the flip angle is 15°, the matrix size is 288 × 224 × 176, the voxel size is 1 mm × 1 mm × 1 mm, and the receive bandwidth is 250 Hz / pixel. The acceleration factor is 2, and the scanning time is 7 minutes and 44 seconds.

[0032] Step 2: Preprocess the modulus diagram and the phase diagram of all echoes to obtain the region of interest and the total field image;

[0033] Step 3: Calculate the surface triangular model according to the region of interest obtained in step 2;

[0034] Step 4: Use the surface triangular model obtained in step 3 to construct the optimization formula of the total field inversion in combination with other sequence-related parameters, which is expressed as the following equation:

[0035] (1)

[0036] Wherein, is the magnetization value to be solved, is the boundary condition to be solved, for intermediate variables, For preprocessing, For magnitude weighting, It is a dipole nucleus. For the integral operator of the Laplace equation, The total field obtained in step 2, It is the regularization weight. It is a structural information mask obtained from the model diagram.

[0037] In formula (1), the operator The following system of integral equations, constructed using Green's function, is obtained:

[0038] (2)

[0039] in, is the Green's function of the Laplace equation, S is the brain surface, and n′ is the unit outward normal vector of S'. In this method, the algorithm uses the brain surface triangulation obtained from step 3 to discretize equation (2) to obtain the operator. Discrete matrix form Used for subsequent numerical calculations. It has the following matrix form:

[0040] (3)

[0041] Formula (3) is The specific composition consists of four matrices: T, Q, T', and Q', with elements being the integrals of the surface triangular model faces. -1 The expression for each element of T, Q, T', Q' is also given in formula (3). The subscripts i and j represent the numbers of the surface triangles, and (x, y, z) are the coordinates in the three-dimensional image space, used to represent the numbers of the matrix elements.

[0042] Step 5: Use the fast multipole expansion algorithm to store the matrix. The 3D field map matrix and the surface triangular mesh are organized into an octree structure. The center of the octree is aligned with the center of the field map matrix to ensure that the center position of each level of cube is consistent with the matrix coordinate system. The cubes are recursively subdivided starting from the root node covering the entire domain; the voxels / sampling points of the field map matrix are assigned to the corresponding cubes according to their spatial coordinates; the boundary triangles are divided into the corresponding octree sub-cubes according to their centroid positions, and the center coordinates, side length / radius, and level of each node are calculated accordingly. Then, the Green's function is expanded spherically harmonicly with the nodes as the center as follows:

[0043] (4)

[0044] where is the expansion node, r represents the expansion start coordinate, and r' represents the expansion end coordinate. and are functions related to the expansion start and end coordinates, and specific forms are given below. is the associated Legendre polynomial. and are the orders of the Legendre polynomials. and are spherical coordinate parameters of the corresponding vector according to the expansion. In each subsequent step of L calculation, an uplink transmission from leaf to root and a downlink transmission from root to leaf along the octree are required to perform the far-field calculation, and a direct matrix multiplication calculation is performed on the near-field pairs determined by the tree's neighbor list. Finally, the far-field and near-field contributions are summed.

[0045] Step 6: Numerically calculate the solution of formula (1) in step 4 using the Gauss-Newton method to obtain the magnetic susceptibility distribution image. The specific loop calculation implementation is as follows:

[0046] 1) Set the to-be-solved magnetic susceptibility related intermediate variables and initial values and ;

[0047] 2) For the nth step, use the discrete matrix L instead of operator formula (1) at and to perform the following Taylor expansion: (5)

[0048] where is the loss function of the nth step;

[0049] 3) Find when the gradient is 0 and , that is, solve the following equation:

[0050] (6)

[0051] where is a parameter set to prevent the denominator from being 0;

[0052] 4) Calculate and , and repeat steps 2) and 3) until and are less than the threshold value to obtain the final result.

[0053] The magnetic susceptibility distribution map in this example can be obtained according to the above steps, as shown in Figure 3 FIG. 6.

Claims

1. A method for calculating a quantitative magnetic susceptibility image based on total field inversion using Green's function, characterized in that, The method includes the following steps: Step 1: Acquire raw magnetic resonance images using a multi-echo gradient echo sequence; Step 2: Preprocess the mode map and phase map of all echoes to obtain the region of interest and total field image; Step 3: Calculate the surface triangulation model based on the region of interest obtained in Step 2; Step 4: Using the surface triangular model obtained in Step 3 combined with the modal information obtained from the magnetic resonance sequence in Step 1, construct the optimized formula for total field inversion, expressed as the following equation: （1）； in, The desired magnetic susceptibility value is... The boundary conditions are to be determined. for intermediate variables, For preprocessing, For magnitude weighting, It is a dipole nucleus. For the integral operator of the Laplace equation, The total field obtained in step 2, It is the regularization weight. It is a structural information mask obtained from the model diagram; In formula (1), the operator The following system of integral equations, constructed using Green's function, is obtained: （2）； in, It is the Green's function of the Laplace equation, S is the brain surface, and n′ is the unit outward normal vector of S'; the operator is obtained by discretizing formula (2) using the brain surface triangulation model obtained from step 3. Discrete matrix form Used for subsequent numerical calculations. It has the following matrix form: （3）； Formula (3) is The specific composition consists of four matrices: T, Q, T', and Q', with elements being the integrals of the surface triangular model faces. -1 The expression for each element of T, Q, T', Q' is also given in formula (3); where the subscripts i and j are the numbers of the surface triangles, and (x,y,z) are the coordinates in the three-dimensional image space, used to represent the number of the matrix elements; Step 5: Compress the matrix using the fast multipole expansion algorithm. The size is small enough to be stored on a desktop workstation or a small server. The algorithm first constructs an octree mesh for the 3D image matrix and surface triangular model and calculates the mesh parameters; then the Green's function is spherically harmonicly expanded with the mesh nodes as the center. Step 6: Use the Gauss-Newton method to numerically calculate the solution of formula (1) in step 4 and obtain the magnetic susceptibility distribution image.

Images