A method for eliminating large area water-fat interchanges in magnetic resonance water-fat separation techniques

Abstract

This invention discloses a method for eliminating large-area water-lipid exchange in magnetic resonance imaging (MRI) water-lipid separation technology, relating to the field of medical image processing technology. The method includes: inputting original images with multiple echoes and obtaining downsampled images using different downsampling factors; constructing three initial field maps and performing dimensionality reduction simplification on each of the three initial field maps; then iterating the field maps based on graph cuts to obtain the corresponding field maps and their ranges; continuously comparing and judging the ranges of the three field maps to determine the final field map; finally, substituting the final field map into the signal intensity model formula to solve for the water-lipid separation image. This method can obtain clear and accurate water-lipid separation images of various parts of the human body, reducing the occurrence of large-area water-lipid exchange; suppressing adipose tissue signals and increasing image contrast; enhancing the scanning effect; and obtaining fat fractions to reflect the development of osteoporosis and fatty liver.

Description

Technical Field

[0001] This invention relates to the field of medical image processing technology, specifically to a method for eliminating large-area water-lipid exchange in magnetic resonance water-lipid separation technology. Background Technology

[0002] Currently, the magnetic resonance fat suppression techniques used in clinical practice include: frequency-selective saturation (FS), short inversion time inversion recovery (STIR), frequency-selective inversion pulse fat suppression techniques (including SPIR and SPAIR techniques), selective water or fat excitation techniques (including PROSET, WATS, and SPGR), and chemical shift inversion imaging techniques (Dixon).

[0003] FS (Fat Suppression) utilizes the difference in precession frequencies of protons in fat and water molecules, applying a prepulse with the same frequency as fat protons to saturate the protons in fat, thus achieving fat suppression. Its advantages include high selectivity, suppressing almost entirely adipose tissue signals with minimal impact on signals from other tissues. However, it requires high field strength and magnetic field uniformity, and is therefore generally used with medium- to high-field machines. STIR (Speed-Inverted Radiation) is an amplitude-selective saturation technique that images based on different T1 values ​​of human tissues. It has low field strength dependence and field uniformity requirements, but poor image signal-to-noise ratio. Frequency-selective inversion pulse (SPIR) fat suppression technology can be seen as a combination of FS and STIR. SPIR shortens the scan time but requires higher field strength and field uniformity, while SPIR provides more thorough fat suppression but prolongs the scan time. PROSET and WATS technologies utilize the chemical shift effect between fat and water, possessing spatial selectivity in addition to frequency selection, selectively exciting either fat or water. SPGR technology eliminates T2 interference by applying a scrambling RF pulse, thus reflecting T1 better.

[0004] Dixon first proposed a chemical shift MRI method for separating water and fat in 1984. The Dixon technique utilizes the difference in resonance frequencies between water and fat to obtain in-phase, out-of-phase, water, and fat maps in a single pass, reducing the field strength requirement but making it susceptible to motion artifacts. However, improved Dixon techniques (FTED, IDEAL) significantly shorten scan time, reduce respiratory motion artifacts, and completely separate water and fat, showing broad application prospects. In 1997, Xiang Qingsan proposed a new three-point Dixon method with direct phase encoding, sampling three points (0, π / 2, π) to determine which image is water and which is fat while separating them. In 2007, Berglund and Kullberg showed that by considering only two periodically repeating candidates for the B0 partial resonance of each voxel, the water-fat separation problem can be solved non-iteratively using a single so-called quadratic pseudo-Boolean optimization (QPBO) map cut. In 2010, Hernando proposed a water-fat separation method based on a graph cut algorithm. This method first discretizes the optimization problem, then sets iterative conditions, and solves subproblems using graph cut until a global minimum is found. However, this algorithm still yields incorrect water-fat separation results in low signal-to-noise ratio regions. In 2016, Cheng Chuanli introduced phasor estimation into the solution model, avoiding phase untangling of the field map. His proposed seed pixel identification method and region growing method are executed independently at different resolutions. Multiple phasor maps are obtained at lower resolutions and then merged into a new seed map, which is used to generate the final phasor map at the optimal resolution. Based on this phasor map, the final fat and water images are reconstructed. The self-feedback mechanism of phasor estimation ensures the reliability of seed pixel selection at the optimal resolution.

[0005] However, current graph cut algorithms are prone to producing incorrect water-lipid separation results when adipose tissue is not continuous, especially in complex areas of the human body such as the pelvis and cervical spine. Therefore, this invention proposes a method for eliminating large-area water-lipid exchange in magnetic resonance water-lipid separation technology. Summary of the Invention

[0006] The purpose of this invention is to provide a method for eliminating large-area water-fat exchange in magnetic resonance imaging (MRI) water-fat separation technology. By combining a fat multi-peak model and variable dimensionality reduction (and also utilizing compressed sensing, parallel imaging, and half-Fourier acquisition technology), it can greatly accelerate the scanning sequence, suppress motion artifacts caused by physiological factors and patient movement, reduce scanning time, and bring a better examination experience to doctors and patients. This represents the future development direction of the field of magnetic resonance imaging.

[0007] The technical problem solved by this invention is: how to separate the fat and water components in medical images to avoid large-scale water-fat exchange.

[0008] This invention can be achieved through the following technical solution: a method for eliminating large-area water-lipid exchange in magnetic resonance water-lipid separation technology, the method comprising the following steps:

[0009] Step 1: Input the original images of N echoes, where N≥3, and use different downsampling factors to obtain downsampled images;

[0010] Step 2: Construct three initial field maps, labeled as field map fm1Initial, fm2Initial, and field map fm3Initial respectively. Field maps fm1Initial and fm2Initial are obtained by taking the conjugate subtraction result of the image domain data of their first two echoes. Field map fm3Initial is the result after setting the initial value A on the blank matrix points.

[0011] Step 3: For each of the three initial field maps, reduce their dimensionality and simplify them to binary to select a model. Then, iterate the field maps based on graph cuts to obtain the corresponding field maps fm1, fm2 and fm3, and calculate the range of the three obtained field maps.

[0012] Step 4: Perform continuous comparison and judgment on the ranges of the three field maps, and determine the final field map fm;

[0013] Step 5: Substitute the final field image fm into the formula for the single-point q signal intensity model of the field at different times, and then solve for the water-fat separation image.

[0014] A further technical improvement of this invention lies in: during the chemical shift-based sequence scanning process, different TE times (t1, t2, ..., t...) are collected. n The signal for t n At time q, the signal strength s of a single point q q (t n ) is represented as:

[0015]

[0016] Where: f B,q (Unit: Hz) is due to local magnetic field shift caused by non-uniformity of the static magnetic field; ρ W,q and ρ F,q These represent the amplitude values ​​for water and fat components, respectively.

[0017] f F (Unit: Hz) represents the frequency shift of fat relative to water, and is a known parameter;

[0018] and This represents the time constant of the transverse magnetization vector decay.

[0019] A further technical improvement of the present invention lies in the following specific steps for dimensionality reduction and simplification of the initial field diagram:

[0020] (1) Construct an error function for point q and minimize it to obtain the maximum likelihood estimates of the water or fat amplitude ρ and the local magnetic field offset f:

[0021]

[0022] s q =[s q (t1)···s q (t n )]

[0023] (2) Using the penalized maximum likelihood framework, calculate the estimate of the complete field map:

[0024]

[0025] Where, δ q It is a local nearest neighbor of point q, μ is a regularization coefficient that balances data consistency and solution smoothness, and W q,j It is a spatially dependent weight, V(f) B,q ,f B,j ) is the penalty term for field map roughness, and V(f) B,q ,f B,j )=(f B,q ,f B,j ) 2 ;

[0026] (3) Dimensionality reduction is achieved using variable projection, by... Minimize to obtain the nonlinear parameter f B,q , where ψ(f B,q ) is an N×2 matrix, and n = 1, 2, ..., N + Represents the inverse pseudomatrix;

[0027] Therefore, the field map estimation in step (2) can be discretized as follows:

[0028]

[0029] (4) By limiting the discrete range, there are Ω possibilities for the field at each point. Taking a subset of these possibilities, such that the field map at each point has only two possibilities—the current field map value and the field map value that will change—the field map expression in step (3) can be simplified to:

[0030]

[0031] A further technical improvement of the present invention lies in: using the graph cut algorithm for the simplified field graph expression, making pairwise iterative selections between local minima, determining whether to use the current field value or the next field value, and finding the minimum value within the global scope, so as to determine the field graphs fm1, fm2, and fm3 corresponding to the three initial field graphs and the corresponding field ranges range1, range2, and range3.

[0032] A further technical improvement of the present invention lies in: the judgment process for determining the final field graph is as follows:

[0033] When range1 ≤ range2 and range1 ≤ range3, let the final field graph fm = fm1, otherwise execute the next judgment;

[0034] When range2 ≤ range1 and range2 ≤ range3, let the final field graph fm = fm2, otherwise execute the next judgment;

[0035] Judge whether |range1 - range3| < M is satisfied. If so, let the final field graph fm = fm1;

[0036] If not, let the final field graph fm = fm3, where M is the difference in the field range obtained through experiments.

[0037] A further technical improvement of the present invention lies in: for the next field value, it can be selected from the following three different sets:

[0038] Γ β :

[0039] Γ + :

[0040] Γ - :

[0041] Among them, is the next possible field value of point q, β is a constant, is a series of local minima of R0(f B ; s q ) at point q; Γ β means that each point on the field graph makes a uniform jump, and Γ + and Γ - represent jumps corresponding to voxels.

[0042] A further technical improvement of the present invention lies in: the initial value A is specifically set to 150 Hz.

[0043] Compared with the prior art, the present invention has the following beneficial effects:

[0044] This method can obtain clear and accurate water-lipid separation images of various parts of the human body (head, spine, breast, liver, and limb joints), reduce the occurrence of large-area water-lipid exchange, suppress adipose tissue signals, increase image contrast, enhance scanning effects, and obtain fat fractions to reflect the development of osteoporosis and fatty liver. Attached Figure Description

[0045] To facilitate understanding by those skilled in the art, the present invention will be further described below with reference to the accompanying drawings.

[0046] Figure 1 This is a flowchart illustrating the method execution of the present invention.

[0047] Figure 2 This is a schematic diagram illustrating the possible selection of field values ​​for the corresponding point in this invention.

[0048] Figure 3 This is a schematic diagram of the graph cut problem logic of the present invention;

[0049] Figure 4 This is a comparison chart of the water-lipid separation results obtained from three different field diagrams of the present invention;

[0050] Figure 5 This is a comparison chart of the results of the original graph cut algorithm and the method of the present invention in the sagittal plane image of the pelvis. Detailed Implementation

[0051] To further illustrate the technical means and effects of the present invention in achieving its intended purpose, the following detailed description of the specific implementation methods, structures, features, and effects of the present invention, in conjunction with the accompanying drawings and preferred embodiments, is provided.

[0052] Please see Figure 1 As shown, a method for eliminating large-area water-lipid exchange using magnetic resonance water-lipid separation technology includes the following steps:

[0053] Step 1: First, input the original image of N echoes, where N≥3. Use different downsampling rates to obtain the downsampled image. The downsampling factors used in this invention are Q and 2Q, that is, Q times and 2Q times the sampling, and the downsampling factor Q is 2, 3, 4, etc.

[0054] Step 2: Take the image domain data of the first two echoes of each and perform conjugate subtraction to obtain the initial field maps fm1Initial and fm2Initial. At the same time, create a blank matrix and set the initial value of each point in the matrix to A as the field map fm3Initial. The initial value A of the points in the field map is set to the value within the initial estimated field range. For example, if the initial estimated field range is [-300Hz, 300Hz], then A = 150Hz can be set.

[0055] Step 3: Reduce the dimensionality of each of the three initial field map values and simplify them to a binary selection model. Subsequently, perform iterative field map based on graph cuts to obtain the corresponding field maps fm1, fm2, and fm3, and calculate the ranges of the three obtained field maps, which are range1 = max(fm1) - min(fm1), range2 = max(fm2) - min(fm2), and range3 = max(fm3) - min(fm3).

[0056] Step 4: Judge the field map ranges in Step 3:

[0057] If range1 ≤ range2 and range1 ≤ range3, then let the final field map fm = fm1; otherwise, perform the next judgment;

[0058] If range2 ≤ range1 and range2 ≤ range3, then let the final field map fm = fm2; otherwise, continue to judge whether |range1 - range3| < M is satisfied: if yes, then let the final field map fm = fm1; if not, then let the final field map fm = fm3;

[0059] Among them, the field range difference M is a frequency value obtained through experimental testing and can be set to 50 Hz;

[0060] Step 5: Substitute the final result of the field map into the single-point q signal intensity model formula of the corresponding field at different times, and then solve the water-fat separation image.

[0061] Specifically, during the sequence scanning based on chemical shift, signals at different TE times (t1, t2, ···, t n ) are collected. For the t[[ID=2​​​​​​​​​​​​​​​​​​​​​​​​After magnetic resonance excitation, the longitudinal magnetization vector increases exponentially during relaxation, with a time constant of T1, while the transverse magnetization vector decreases exponentially during relaxation, with a time constant of T2. In reality, the main magnetic field B0 cannot be perfectly uniform because the rotation frequency of hydrogen atoms is related to the strength of B0. Non-uniform B0 leads to different rotation frequencies of hydrogen atoms at different locations, resulting in asynchronous rotation. This accelerates the decay of the transverse magnetization vector, which also occurs exponentially, with a time constant of T2. and It is smaller than T2.

[0067] Minimizing the following error function yields the maximum likelihood estimates of ρ and f:

[0068]

[0069] s q =[s q (t1)···s q (t n (2)

[0070] By invoking the penalized maximum likelihood framework, we can compute an estimate of the complete field diagram:

[0071]

[0072] Where, δ q It is a local nearest neighbor of point q, μ is a regularization coefficient that balances data consistency and solution smoothness, and W q,j It is a spatially dependent weight, V(f) B,q ,f B,j ) is the penalty term for field map roughness, and V(f) B,q ,f B,j )=(f B,q ,f B,j ) 2 .

[0073] By using variable projection (VARPRO) to reduce the dimensionality of the problem, the nonlinear parameter f B,q It can be obtained by minimizing the following formula:

[0074] Wherein, ψ(f B,q ) is an N×2 matrix, and n = 1, 2, ..., N + This represents the inverse pseudomatrix.

[0075] Therefore, the field map estimation can be discretized as follows:

[0076]

[0077] The field diagram can be discretized to a range of ±1500Hz, with each spacing of 2–4Hz. Therefore, the field at each point has Ω possible values. When the spacing is 3Hz, Ω = 1500 / 3 = 500. The above equation then becomes:

[0078]

[0079] Taking a subset Γ of Ω, the field pattern at each point has only two possibilities: the current field pattern value and the field pattern value to be changed. The above formula can be further simplified to:

[0080]

[0081] The field graph expression can be solved quickly using the graph cut algorithm. Initially, the field of all points is 0. Then, in continuous iteration, it is determined whether to use the current field or the next field. The optimal solution that satisfies the above expression can be found globally.

[0082] Using three different sets Γ, the optimal solution can be found quickly:

[0083] Γ β :

[0084] Γ + :

[0085] Γ - :

[0086] in, Let q be the next possible field value, and β be a constant. It is a series of R0(f) B ;s q The local minimum at point q; Γ β This indicates that each point on the field diagram must jump uniformly, Γ + and Γ - This indicates a jump corresponding to voxel-related behavior.

[0087] like Figure 2 As shown, R0(f) represents a single point. B ;s q There are many local minima, but by iterating and jumping around, we can find the final global minimum.

[0088] like Figure 3 The diagram shown illustrates a graph cut problem, where the points are v1, v2, ..., v Q s represents the source node, t represents the sink node, and d represents the sink node. 1t ,d12 ,d s1 ,...,d qt This represents the weight between points;

[0089] It should be noted that, for the points and weights in the graph cut, the field graph of each point is discretized into L possible values. Therefore, in a K x K graph, the field graph has L possible values. KxK For example, in an image with 192x192 pixels, if each pixel has 500 possible field map values, then the total number of possible field map values ​​is 500. 192*192 =500 36864 In a graph cut, there are a total of K+2 points, with the two extra points being the source and sink. The weights of the lines connecting the points consist of regularization weights and residual weights.

[0090] In a network flow, the maximum flow from the source to the sink is equal to the minimum capacity of the set of edges that would cause the network flow to be interrupted if removed from the network. Traverse all possible paths, and use the minimum edge in each path as the maximum flow of that path. Subtract this flow to obtain the edges with no flow, which are then considered cuts. This algorithm can solve the dimensionality-reduced optimization problem, thus obtaining the true field graph. Finally, by substituting into formula (1), the image of water-lipid separation can be obtained.

[0091] like Figure 4 As shown, the water-lipid separation results obtained from three different field maps are compared. The results of field map fm1 are used to avoid large-area water-lipid exchange.

[0092] like Figure 5 The figure shows a comparison of the results of the original graph-cut algorithm and the algorithm proposed in this invention in the sagittal plane image of the pelvis. The original graph-cut algorithm performs poorly in the complex pelvic structure, with only 6 out of 18 layers being correctly separated (Figures a and b are the lipid and water maps obtained by the original graph-cut algorithm, respectively). In contrast, the algorithm proposed in this invention achieves good lipid and water separation results in both the left and right planes (Figures c and d are the lipid and water maps obtained by the algorithm proposed in this invention, respectively).

[0093] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some modifications or alterations to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent changes and alterations made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the scope of the present invention.

Claims

1. A method for eliminating large area water-fat interchanges in magnetic resonance water-fat separation techniques, characterized by, The method includes the following steps: Step 1: Input the original images of N echoes, where N≥3, and obtain the downsampled images with different downsampling factors; Step 2: Construct three initial field maps, respectively labeled as field map fm1 Initial, fm2 Initial, and field map fm3 Initial. Among them, for field map fm1 Initial and field map fm2 Initial, take the conjugate subtraction results of the image domain data of their respective first two echoes, and field map fm3 Initial is the result after setting the initial value A at the blank matrix points; Step 3: Perform dimensionality reduction and simplify each of the above three initial field maps to a binary selection model, and then perform field map iteration based on graph cuts to obtain the corresponding field maps fm1, fm2, and fm3, and calculate the ranges range1, range2, and range3 of the three field maps; Step 4: Continuously compare and judge the ranges of the three field maps, and determine the final field map fm; Step 5: Substitute the final field map fm into the single-point q signal intensity model formula of the corresponding field at different times, and then solve the water-fat separation image; The judgment process for determining the final field map is as follows: When range1≤range2 and range1≤range3, let the final field map fm = fm1, otherwise perform the next judgment; When range2≤range1 and range2≤range3, let the final field map fm = fm2, otherwise perform the next judgment; Judge whether |range1 - range3| < M is satisfied. If so, let the final field map fm = fm1; If not, let the final field map fm = fm3, where M is the range difference of the field obtained from the experiment.

2. The method for eliminating large-area water-fat exchange in magnetic resonance water-fat separation technology according to claim 1, characterized in that, In a chemical shift based sequence scanning procedure, signals are acquired at different TE times for a time instant, the signal intensity of a single point q is expressed as: ； in: The unit is Hz, which represents the local magnetic field shift caused by the inhomogeneity of the static magnetic field; and These represent the amplitude values ​​for water and fat components, respectively. The unit is Hz, which represents the frequency shift of fat relative to water, and is a known parameter; = ,and This represents the time constant of the decay of the transverse magnetization vector of the imaged object at point q.

3. The method for eliminating large-area water-lipid exchange in magnetic resonance water-lipid separation technology according to claim 2, characterized in that, The specific steps for dimensionality reduction and simplification of the initial field map include: (1) Construct an error function formula for point q and perform minimization to obtain the maximum likelihood estimates of the water or fat amplitude ρ and the local magnetic field offset f; (2) Use the penalized maximum likelihood framework to calculate the estimation of the field map; ； in, It is a point Local neighbors, It is a regularization coefficient that balances data consistency and solution smoothness. It is a space-dependent weight. It is a penalty term for field map roughness, and ; (3) Dimensionality reduction is achieved using variable projection, by... Minimize to obtain nonlinear parameters ,in, It is The matrix, and , , , Represents the inverse pseudomatrix; Thus, the field map estimation in step (2) can be discretely expressed as: ； Limit the discrete range, then there are Ω possibilities for the field of each point. Take its subset so that there are only two possibilities for the field map of each point, namely the current field map value and the field map value to be changed. Then the field map expression in step (3) can be simplified as: 。 4. The method for eliminating large-area water-lipid exchange in magnetic resonance water-lipid separation technology according to claim 3, characterized in that, Use the graph cut algorithm for the simplified field map expression, perform pairwise iterative selection among local minima, determine whether to use the current field value or the next field value, and find the minimum value within the global range, so as to determine the field maps fm1, fm2, and fm3 corresponding to the three initial field maps and the corresponding field ranges range1, range2, and range3.

5. A method for eliminating large-area water-lipid exchange in magnetic resonance water-lipid separation technology according to claim 4, characterized in that, For the next field value mentioned above, it can be selected from the following three different sets: ； in, It is the next possible field value of point q. It is a constant. It is a series The local minimum at point q; This means that each point on the field diagram must jump uniformly. and This indicates a jump corresponding to a voxel.

6. The method for eliminating large-area water-lipid exchange in magnetic resonance water-lipid separation technology according to claim 1, characterized in that, The specific value of the initial value A is set to 150Hz.

Images