An optimized method and system for phase unwrapping guided by a quality map

Abstract

The present application relates to an optimized algorithm and system for phase unwrapping guided by a quality map, and relates to the technical field of image processing. It solves the problem that when the phase image to be processed contains strong noise, rapid phase changes, or discontinuous regions, the true phase difference between adjacent pixels may be greater than π. In this case, it is very difficult to accurately recover the true phase from the wrapped phase image. The method includes: calculating a phase derivative variance quality map based on the MRI image; guided by the phase derivative variance quality map, applying a phase unwrapping algorithm to start unwrapping from the point with the lowest phase quality value and transition to high-quality regions until all high-quality regions are unwrapped to obtain a phase map. The present application has the following effects: improving the accuracy of phase unwrapping for MRI images, thus better ensuring the accuracy of the obtained MRI phase map and better meeting the requirements for real-time display of thermograms.

Description

Technical Field

[0001] The present application relates to the field of image processing technology, and in particular, to an optimization method and system for phase unwrapping guided by a quality map. Background Art

[0002] Magnetic resonance imaging (MRI) technology has been widely used in clinical diagnosis due to its advantages such as no ionizing damage caused by radioactivity.

[0003] In terms of clinical diagnosis, magnetic resonance imaging (MRI) has high tissue resolution and has unique advantages in preoperative lesion localization. The radiofrequency ablation technology guided by MRI can perform ablation treatment on lesions such as tumors. The key part in radiofrequency ablation is the control and real-time display of the heating temperature of the lesion tissue. Since the thermogram is proportional to the MRI phase map, it is necessary to perform phase unwrapping on the MRI image to calculate the thermogram in real time.

[0004] It is known that during the process of phase unwrapping of MRI images, due to the periodicity of the phase, the phase data obtained from magnetic resonance imaging is restricted to the principal value interval of (-π, π]. Therefore, it is necessary to decouple the wrapped phase to obtain the correct phase map.

[0005] In related technologies, the main phase unwrapping algorithms for phase are as follows: path tracking; cost function optimization scheme; Markov random field algorithm; least squares method; minimum spanning tree method; using Poisson equation for unwrapping, etc. The above algorithms are usually based on the assumption that the true phase is smooth and the phase difference between any two adjacent points is not greater than π.

[0006] In view of the above related technologies, the inventor found the following defects: The prerequisite for the above unwrapping algorithms to achieve unwrapping is that the true phase difference between adjacent pixels is less than π. In fact, when the phase image to be processed contains strong noise, rapid phase changes, or disconnected regions, the true phase difference between adjacent pixels may be greater than π. In this case, it is very difficult to accurately recover the true phase from the wrapped phase image. Summary of the Invention

[0007] In order to improve the accuracy of phase unwrapping of MRI images, thereby better ensuring the accuracy of the obtained MRI phase map and better meeting the need for real-time display of the thermogram, the present application provides an optimization method and system for phase unwrapping guided by a quality map.

[0008] In the first aspect, the present application provides an optimization method for phase unwrapping guided by a quality map, adopting the following technical solution:

[0009] An optimization method for phase unwrapping guided by a quality map includes:

[0010] The phase derivative variance quality map is calculated based on the MRI image;

[0011] Guided by the phase derivative variance quality map, the phase unwrapping algorithm is applied to start unwrapping from the point with the lowest phase quality value and transition to the high-quality region until all high-quality regions are unwrapped, obtaining a phase map.

[0012] By adopting the above technical solution, since the phase derivative variance quality map is selected as the phase guidance map for guidance, which utilizes the phase gradient variance information within the region, it can overcome the drawback that the conventionally selected maps such as the pseudo-correlation quality map and the maximum phase gradient quality map are prone to misidentifying reliable regions with large phase gradients but no undersampling and no noise as low-quality regions. Under the guidance of the phase derivative variance quality map, and considering that the phase transformation of the tissue part in the MRI image is relatively gentle and the phase difference in the background part is relatively large, the method of starting unwrapping from the low phase quality number and extending to the high-quality region can complete the unwrapping of all regions with a relatively high unwrapping accuracy, obtaining a relatively accurate phase map.

[0013] Optionally, calculating the phase derivative variance quality map based on the MRI image includes:

[0014] An initial phase map is calculated based on the MRI image. Here, the MRI image is defined as IMG, with height M and width N. IMG (m，n) represents taking the value at the m-th row and n-th column of the matrix IMG, where both m and n are integers; the neighborhood set of the point with coordinates (m, n) and its upper, lower, left, and right coordinates is set as follows:

[0015] C{(m, n)} = {(m, n), (m + 1, n), (m - 1, n), (m, n - 1), (m, n + 1)};

[0016] Calculate the phase gradient map in the X direction;

[0017] Calculate the phase gradient map in the y direction;

[0018] Calculate the neighborhood mean of the (m, n) point in the X direction phase gradient map;

[0019] Calculate the neighborhood mean of the (m, n) point in the y direction phase gradient map;

[0020] Calculate the phase quality map Q, with height M and width N.

[0021] By adopting the above technical solution, the specific formula algorithm for analyzing the phase derivative variance quality map based on the MRI image is specifically disclosed. The parallel processing of some parallel step algorithms effectively shortens the calculation of the phase derivative variance quality map.

[0022] Optionally, the initial phase map is calculated based on the MRI image as follows: P (m，n) = arctan(IMG (m，n) ) where m ∈ [0, M), n ∈ [0, N), and P (m，n) is the initial phase map.

[0023] By adopting the above technical solution, the wrapped phase values within one period phase interval can be extracted from the MRI image using the arctangent function, thereby forming an initial phase map that meets the requirement of the principal value interval of the phase value being restricted to (-π, π].

[0024] Optionally, the calculation of the neighborhood mean value of the (m, n) point on the X - direction phase gradient map is as follows:

[0025]

[0026] By adopting the above technical solution, the calculation of the neighborhood mean value of the (m, n) point on the X - direction phase gradient map is specifically disclosed.

[0027] Optionally, the calculation of the neighborhood mean value of the (m, n) point on the y - direction phase gradient map is as follows:

[0028]

[0029] By adopting the above technical solution, the calculation of the neighborhood mean value of the (m, n) point on the y - direction phase gradient map is specifically disclosed.

[0030] Optionally, the calculation of the phase quality map Q includes:

[0031] Initialized to the maximum value: Q (m，n) = 100 where m ∈ [0, M), n ∈ [0, N);

[0032]

[0033] By adopting the above technical solution, the variances of the differences in the horizontal and vertical directions are respectively statistically calculated and then added together. It represents the deviation range of the difference at a certain point from its mathematical expectation. Therefore, it is very sensitive to the discontinuity of the phase and at the same time has no noise problem of the correlation coefficient map.

[0034] Optionally, the phase unwrapping algorithm includes:

[0035] Initialization: Let the unwrapped phase map be a matrix IM_uwp of M * N, the neighborhood marking map be a matrix Map_adj of M * N, the unwrapped marking map be a matrix Map_uwp of M * N, and the four - neighborhood reference point set be G, specifically including the following:

[0036] Calculate the coordinates of the minimum - value point, Initialize the neighborhood marking graph to 0, and mark the neighborhood of point (mm, nn) as 1:

[0037] Map_adj (mm-1，nn) = 1; Map_adj (mm+1,nn) = 1: Map_adj (mm，nn-1) = 1; Map_adj (mm，nn+1) = 1;

[0038] Initialize the value of the unwrapped point graph to 0, and mark point (mm, nn) as 1: Map_uwp (mm，nn) = 1; Initialize the value of the phase graph after unwrapping to 0, and calculate the phase value of point (mm, nn): IM_uwp (mm，nn) = P (mm，nn) ; Initialize set G to be an empty set; and,

[0039] Update the quality graph Q, calculate the minimum coordinate, update set G, and select the neighborhood calculation mode mod;

[0040] Calculate the phase graph IM_uwp after unwrapping according to the calculation mode classification, and update the neighborhood graph Map_adj and the marking graph Map_uwp;

[0041] Repeat the operations of updating the quality graph Q, calculating the minimum coordinate, updating set G, and selecting the neighborhood calculation mode mod, until the neighborhood marking graph Map_adj is all 0, then all points are unwrapped.

[0042] By adopting the above technical solution, considering that the phase transformation of the tissue part in the magnetic resonance image is relatively gentle, and the phase difference of the background part is relatively large, start unwrapping from the selected phase point with the lowest quality, then unwrap the nearby pixel points, put the unwrapping result into the adjacent array, select the point with the lowest quality from the above array, and delete it, then connect this point and put its nearby pixel points into the nearby table; finally, sort the points in the nearby table in ascending order of quality, and repeat the operation continuously to continuously expand the phase unwrapping area, and finally calculate the point with the lowest quality and the highest quality.

[0043] Optionally, the phase unwrapping algorithm includes:

[0044] Initialization, define the data structure: Let the phase graph after unwrapping be a matrix IM_uwp of M*N, the neighborhood data structure be a red-black tree RBTree, the unwrapped point marking graph be a matrix Map_uwp of M*N, and the set of four-neighborhood reference points be G. Specifically as follows: Calculate the coordinates of the minimum point, Initialize the neighborhood tree, insert the neighborhood nodes of the minimum point (mm, nn), RBTree→insert({Q(mm - 1，nn) , (mm - 1, nn)});

[0045] RBTree->insert({Q (mn+1,nn) , (mm + 1, nn)});

[0046] RBTree->insert({Q (m，nn-1) , (mm, nn - 1)));

[0047] RBTree->insert({Q (m，nn+1) , (mm, nn + 1)}); Initialize the unwrapped marker graph to 0 and mark the point (mm, nn) as 1: Map_uwp (mm，nn) = 1; Initialize the value of the unwrapped phase graph to 0 and calculate the phase value of the point (mm, nn): IM_uwp (mm，nn ) = P (mm，nn) ; Initialize the set G to be an empty set; and,

[0048] Update the neighborhood tree RBTree, calculate the minimum coordinate (mi, ni), update the set G, and select the neighborhood calculation mode mod;

[0049] Calculate the unwrapped phase graph IM_uwp according to the calculation mode classification, and update the neighborhood marker tree RBTree and the unwrapped marker graph Map_uwp;

[0050] Repeat the operations of updating the neighborhood tree RBTree, calculating the minimum coordinate (mi, ni), updating the set G, selecting the neighborhood calculation mode mod, and calculating the unwrapped phase graph IM_uwp according to the calculation mode classification, and updating the neighborhood marker tree RBTree and the unwrapped marker graph Map_uwp until the neighborhood tree RBTree is empty, then all points are unwrapped.

[0051] By adopting the above technical solution, which integrates the red - black tree architecture and the global phase unwrapping algorithm, during the process of unwrapping from low - quality phases to high - quality regions, the lowest - quality phase can be found faster, and the overall unwrapping efficiency can be improved on the premise of ensuring the phase unwrapping accuracy.

[0052] Optionally, updating the neighborhood tree RBTree, calculating the minimum coordinate (mi, ni), updating the set G, and selecting the neighborhood calculation mode mod includes:

[0053] Calculate the minimum node of RBTree: {Q (mi，ni) , (mi, ni)} = RBTee->minimum();

[0054] Delete the minimum node of the RBTree: RBTree→erase({Q (mi,ni) , (mi, ni)});

[0055] Update the reference point set G: Calculate the four-neighborhood reference point set G(mi, ni) of the non-boundary point (mi, ni). The upper neighborhood reference point G_up = nozero(Q (mi-1，ni) *Map_uwp (mi-1，ni) );

[0056] The lower neighborhood reference point G_down = nozero(Q (mi+1，ni) *Map_uwp (mi+1，ni) ); The left neighborhood reference point G_left = nozero(Q (mi，ni-1) *Map_uwp (mi，ni-1) ); The right neighborhood reference point G_right = iozero(Q (mi，ni+1 )*Map_uwp (mi，ni+1) ); Then the four-neighborhood reference point set G(mi, ni) = {G_up, G_down, G_left, G_right};

[0057] Select the neighborhood calculation mode: If the point (mi, ni) is a boundary, that is, mi = 0 or mi = M - 1 or ni = 0 or ni = N - 1, then mod = border; if mi is not the upper boundary and the upper neighborhood reference point is the minimum in the four neighborhoods G(mi, ni), that is, mi ≠ 0 and G_up = min(G(mi, ni)), then mod = up; if mi is not the lower boundary and the lower neighborhood reference point is the minimum in the four neighborhoods G(mi, ni), that is, mi ≠ M - 1 and G_down = min(G(mi, ni)), then mod = down; if ni is not the left boundary and the left neighborhood reference point is the minimum in the four neighborhoods G(mi, ni), that is, ni ≠ 0 and G_left = mu(G(mi, ni)), then mod = left; if ni is not the right boundary and the right neighborhood reference point is the minimum in the four neighborhoods G(mi, ni), that is, ni ≠ N - 1 and G_right = min(G(mi, ni)), then mod = right.

[0058] By adopting the above technical solutions, the specific steps of updating the neighborhood tree RBTree, calculating the minimum value coordinates (mi, ni), updating the set G, and selecting the neighborhood calculation mode mod are specifically disclosed.

[0059] In a second aspect, the present application provides an optimized system for phase unwrapping guided by a quality map, adopting the following technical solutions:

[0060] An optimized system for phase unwrapping guided by a quality map, comprising a memory, a processor, and a program stored on the memory and executable on the processor, which can implement an optimized method for phase unwrapping guided by a quality map as described in the first aspect when being loaded and executed by the processor.

[0061] By adopting the above technical solution, through the retrieval of relevant programs, the construction period of the next construction process can be analyzed and obtained more accurately and effectively, and combined with the weather conditions affecting the construction during the construction period of the next construction process, the suitable construction period can be further analyzed and notified to the person in charge of the next construction process in advance, so as to avoid construction in unsuitable weather conditions as much as possible and better ensure the construction quality of the foundation pit support.

[0062] In summary, the beneficial technical effects of this application are as follows: Using the phase derivative variance quality map as the quality guidance map improves the analysis and acquisition efficiency of the quality map due to the parallel computing method. After the quality map guidance, through the optimized algorithm of phase unwrapping, the acquisition efficiency of the MRI phase map is effectively improved, facilitating the timely display of the heat map. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 is the overall flowchart of an optimized method for phase unwrapping guided by a quality map in an embodiment of this application.

[0064] Figure 2 is Figure 1 the specific flowchart of step S100 in.

[0065] Figure 3 is Figure 1 the specific flowchart of one implementation manner of step S200 in.

[0066] Figure 4 is Figure 1 the specific flowchart of another implementation manner of step S200 in.

[0067] Figure 5 is a schematic diagram of an MRI image of the phase to be unwrapped.

[0068] Figure 6 is the phase map of the MRI image during unwrapping.

[0069] Figure 7 is the phase map of the MRI after completion of unwrapping. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0070] The following further describes this application in detail with reference to the accompanying drawings.

[0071] Refer to Figure 1, an optimized method for phase unwrapping guided by a quality map disclosed in this application, includes:

[0072] Step S100, calculating a phase derivative variance quality map based on the MRI image.

[0073] Specifically, the MRI mentioned in step S100 is magnetic resonance imaging. MRI is a phase-sensitive imaging modality. Each pixel value in the complex image after Fourier transform of the MR raw data has a modulus and a phase.

[0074] Among them, the MRI image mentioned in step S100 can be referred to Figure 5 .

[0075] The phase derivative variance quality map mentioned in step S100 is obtained directly from the MRI image by the phase derivative variance method.

[0076] Step S200, guided by the phase derivative variance quality map, applying a phase unwrapping algorithm to start unwrapping from the point with the lowest phase quality value and transition to the high-quality region until all high-quality regions are completely unwrapped to obtain a phase map.

[0077] Specifically, in combination with this application, the phase unwrapping algorithm is explained as follows: The phase data in the complex magnetic resonance array obtained from the MRI image is usually limited to the interval (-π, π]. To obtain the true phase, k 2π values (k is an integer) must be added or subtracted. This process of obtaining the true phase is called phase unwrapping.

[0078] The application of the phase unwrapping algorithm mentioned in step S200 to start unwrapping from the point with the lowest phase quality value and transition to the high-quality region until all high-quality regions are completely unwrapped is further explained as follows: Considering that the phase transformation of the tissue part in the magnetic resonance image is relatively gentle and the phase difference of the background part is relatively large, therefore, the quality of the points in the quality map is sorted in descending order, and the point with the lowest quality in the quality map is unwrapped first, and the point with the highest quality in the quality map is unwrapped last.

[0079] Supplementally, the quality of the above-mentioned points is judged by the height of the pixels. The higher the pixel, the higher the quality of the point.

[0080] Referring to Figure 2 , step S100 includes:

[0081] Step S110, calculating an initial phase map based on the MRI image.

[0082] First, define the MRI image as IMG, with height M, width N, m and n are both integers, IMG (m，n)It represents taking the value of the m-th row and n-th column of the matrix IMG; the neighborhood set of the point with coordinates (m, n) and its adjacent coordinates up, down, left, and right is set as follows: C{(m, n)} = {(m, n), (m + 1, n), (m - 1, n), (m, n - 1), (m, n + 1)}.

[0083] Specifically, the formula for calculating the initial phase map based on the MRI image mentioned in step S110 is as follows: P (m，n) = arctan(IMG (m，n) ) where m ∈ [0, M), n ∈ [0, N), and P (m，n) is the initial phase map.

[0084] Step S1A0, calculate the phase gradient map in the X direction.

[0085] Specifically, the formula for calculating the phase gradient map in the X direction mentioned in step S1A0 is as follows: dx (i，j) = P (i，j+1) - P (i，j) where i ∈ [0, M), j ∈ [0, N - 1).

[0086] Step S1B0, calculate the neighborhood mean of the point (m, n) in the phase gradient map in the X direction.

[0087] Specifically, the calculation of the neighborhood mean of the point (m, n) in the phase gradient map in the X direction mentioned in step S1B0 is as follows:

[0088] Original formula:

[0089] Parallelization: dxu (m，n) = dx (m-1，n) where m ∈ [1, M - 1), n ∈ [1, N - 1);

[0090] dxd (m，n) = dx (m+1，n) where m ∈ [1, M - 1), n ∈ [1, N - 1);

[0091] dxl (m，n) = dx (m，n-1) where m ∈ [1, M - 1), n ∈ [1, N - 1);

[0092] dxr (m，n) = dx (m，n+1) where m ∈ [1, M - 1), n ∈ [1, N - 1);

[0093]

[0094] Step S1a0, calculate the phase gradient map in the y direction.

[0095] Specifically, the formula for calculating the phase gradient map in the y direction mentioned in step S1a0 is as follows:

[0096] dy (i，j) = P (i+1，j) - P (i，j) i ∈ [0, M - 1), j ∈ [0, N).

[0097] Step S1b0: Calculate the neighborhood mean of the point (m, n) in the phase gradient map in the y direction.

[0098] Specifically, the formula for calculating the neighborhood mean of the point (m, n) in the phase gradient map in the y direction mentioned in step S1b0 is as follows:

[0099] Original formula:

[0100] Parallelization: dyu (m，n) = dy (m-1，n) m ∈ [1, M - 1), n ∈ [1, N - 1)

[0101] dyd (m，n) = dy (m+1，n) m ∈ [1, M - 1), n ∈ [1, N - 1)

[0102] dyl (m,n) = dy (m,n-1 )m ∈ [1, M - 1), n ∈ [1, N - 1)

[0103] dyr (m，n) = dy (m，n+1) m ∈ [1, M - 1), n ∈ [1, N - 1)

[0104]

[0105] Step S120: Calculate the phase quality map Q, with height M and width N.

[0106] Specifically, the specific formula for calculating the phase quality map Q, with height M and width N, mentioned in step S120 is as follows:

[0107] Original formula:

[0108] Initialized to the maximum value: Q (m，n) = 100m ∈ [0, M), n ∈ [0, N);

[0109]

[0110] Parallelization:

[0111] dxu (m，n) = dx (m-1，n)m ∈ [1, M - 1), n ∈ [1, N - 1)

[0112] dxd (m，n) = dx (m+1，n) m ∈ [1, M - 1), n ∈ [1, N - 1)

[0113] dxl (m，n) = dx (m，n-1) m ∈ [1, M - 1), n ∈ [1, N - 1)

[0114] dxr (m，n) = dx (m,n+1) m ∈ [1, M - 1), n ∈ [1, N - 1)

[0115]

[0116] Refer to Figure 3 , one implementation of the phase unwrapping algorithm mentioned in step S200 is as follows:

[0117] Step S210, initialization.

[0118] First, ① Define the neighborhood calculation mode: mod = {up, down, left, right, border}

[0119] ② Define the two-point unwrapping function Unwrap(p1, p2):

[0120] Let p1 and p2 be the two points to be unwrapped;

[0121] fix: Round down to zero, round: Round to the nearest integer;

[0122] dp = p2 - p1; dp_corr = dp ÷ 2π;

[0123]

[0124] ③ Define the non-zero initialization function nozero(x):

[0125]

[0126] Secondly, let the unwrapped phase diagram be a matrix IM_uwp of M * N, the neighborhood marking diagram be a matrix Map_adj of M * N, the unwrapped marking diagram be a matrix Map_uwp of M * N, and the four-neighborhood reference point set be G.

[0127] Specifically, step S210 can be further divided as follows: 1. Calculate the coordinates of the minimum value point, 2. Initialize the neighborhood marking diagram to 0, and mark the neighborhood of the point (mm, nn) as 1:

[0128] Map_adj (mm-1，nn) = 1; Map_adj (mm+1，nn) = l; Map_adj (mm，nn-1) = 1; Map_adj (mm，nn+1) = 1; 3. Initialize the value of the unwrapped point map to 0 and mark the point (mm, nn) as 1: Map_uwp (mm，nn) = 1; 4. Initialize the value of the phase map after unwrapping to 0 and calculate the phase value of the point (mm, nn): IM_uwp (mm，nn) = P (mm,nn) ; 5. Initialize the set G to be an empty set.

[0129] Step S220, update the quality map Q, calculate the minimum coordinate, update the set G, and select the neighborhood calculation mode mod.

[0130] Specifically, step S220 can be divided as follows: 1. Update the quality map: Q (m，n) = Q (m，n) * Map_adj (m，n) + 100 * (1 - Map_adj (m，n) ); 2. Calculate the minimum point coordinate after updating the quality map: 3. Update the reference point set G:

[0131] Calculate the four-neighborhood reference point set G(mi, ni) of the non-boundary point (mi, ni); Let the upper neighborhood reference point G_up = nozero(Q (mi-1，ni) * Map_uwp (mi-1，ni) ); The lower neighborhood reference point G_down = nozero(Q (mi+1，ni) * Map_uwp (mi+1，ni) ); The left neighborhood reference point G_left = nozero(Q (mi，ni-1) * Map_uwp (mi，ni-1) ); The right neighborhood reference point G_right = nozero(Q (mi，ni+1) * Map_uwp (mi，ni+1) ); Then the four-neighborhood reference point set G(mi, ni) = {G_up, G_down, G_left, G_right}.

[0132] 4. Select the neighborhood calculation mode:

[0133] ① If the point (mi, ni) is a boundary, that is, mi = 0 or mi = M - 1 or ni = 0 or ni = N - 1, then mod = border.

[0134] ② If mi is not the upper boundary and the upper neighborhood reference point is the smallest in the four-neighborhood G(mi, ni),

[0135] That is, if mi ≠ 0 and G_up = min(G(mi, ni)), then mod = up.

[0136] ③ If mi is not the lower boundary and the reference point of the lower neighborhood is the smallest in the four-neighborhood G(mi, ni), that is, mi ≠ M - 1 and G_down = min(G(mi, ni)), then mod = down.

[0137] ④ If ni is not the left boundary and the reference point of the left neighborhood is the smallest in the four-neighborhood G(mi, ni), that is, ni ≠ 0 and G_left = min(G(mi, ni)), then mod = left.

[0138] ⑤ If ni is not the right boundary and the reference point of the right neighborhood is the smallest in the four-neighborhood G(mi, ni), that is, ni ≠ N - 1 and G_right = min(G(mi, ni)), then mod = right.

[0139] Step S230: Calculate and unwrap the phase diagram IM_uwp according to the calculation mode classification, and update the neighborhood graph Map_adj and the marking graph Map_uwp.

[0140] Specifically, step S230 specifically includes the following:

[0141] 1) mod = border

[0142] ① Phase of the unwrapping point (mi, ni): IM_uwp (mi，ni) = 0;

[0143] ② Update the unwrapped marking graph: Map_uwp (mi，ni) = 1;

[0144] ③ Update the neighborhood marking graph: Map_adj (mi，ni) = 0;

[0145] 2) mod = up

[0146] ① Phase of the unwrapping point (mi, ni): IM_uwp (mi，ni) = Unwrap(IM_uwp (mi-1，ni) , P (mi，ni) );

[0147] ② Update the unwrapped marking graph: Map_uwp (mi，ni) = 1;

[0148] ③ Update the neighborhood marking graph:

[0149] Map_adj (mi，ni) = 0;

[0150] Map_adj (mi+1，ni)= 1 if Map_uwwp(mi + 1, ni) = 0;

[0151] Map_adj (mi，ni-1) = 1 if Map_uwp(mi, ni - 1) = 0;

[0152] Map_adj (mi，ni+1) = 1 if Map_uwp(mi, ni + 1) = 0;

[0153] 3) mod = down

[0154] ① Unwrap the phase of the unwrapping point (mi, ni): IM_uwp (mi，ni) = Unwrap(IM_uwp (mi+1，ni) , P (mi，ni) );

[0155] ② Update the unwrapped marker map: Map_uwp (mi，ni) = 1;

[0156] ③ Update the neighborhood marker map:

[0157] Map_adj (mi，ni) = 0;

[0158] Map_adj (mi-1，ni) = 1 if Map_uwwp(mi - 1, ni) = 0;

[0159] Map_adj (mi，ni-1) = 1 if Map_uwp(mi, ni - 1) = 0;

[0160] Map_adj (mi，ni+1) = 1 if Map_uwp(mi, ni + 1) = 0;

[0161] 4) mod = left

[0162] ① Unwrap the phase of the unwrapping point (mi, ni): IM_uwp (mi，ni) = Unwrap(IM_uwp (mi，ni-1) , P (mi，ni) );

[0163] ② Update the unwrapped marker map: Map_uwp (mi，ni) = 1;

[0164] ③ Update the neighborhood marker map:

[0165] Map_adj (mi，ni) = 0,

[0166] Map_adj (mi-1，ni)= 1 if Map_uwp(mi - 1, ni) = 0;

[0167] Map_adj (mi+1，ni) = 1 if Map_uwp(mi + 1, ni) = 0;

[0168] Map_adj (mi，ni+1) = 1 if Map_uwp(mi, ni + 1) = 0;

[0169] 5) mod = right

[0170] ① Unwrap the phase of the unwrapping point (mi, ni): IM_uwp (mi，ni) = Unwrap(IM_uwp (mi,ni+1) , P (mi，ni) );

[0171] ② Update the unwrapped marker map: Map_uwp (mi，ni) = 1;

[0172] ③ Update the neighborhood marker map:

[0173] Map_adj (mi，ni) = 0, Map_adj (mi+1，ni) = 1 if Map_uwp(mi + 1, ni) = 0;

[0174] Map_adj (mi-1，ni) = 1 if Map_uwp(mi - 1, ni) = 0;

[0175] Map_adj (mi，ni-1) = 1 if Map_uwp(mi, ni - 1) = 0.

[0176] Step S240, repeat Step S220 and Step S230 until the neighborhood marker map Map_adj is all 0, then all points are unwrapped.

[0177] Specifically, the schematic diagram when all points mentioned in Step S240 are unwrapped can be referred to Figure 7 , and the schematic diagram before the unwrapping is completed can be referred to Figure 6 .

[0178] Refer to Figure 4 , another implementation of the phase unwrapping algorithm mentioned in Step S200 is as follows: Step S2A0, initialization.

[0179] First, define the data structure: Let the unwrapped phase diagram be a matrix IM_uwp of M * N, the neighborhood data structure be a red - black tree RBTree, the unwrapped point marker map be a matrix Map_uwp of M * N, and the set of four - neighborhood reference points be G.

[0180] Specifically, step S2A0 includes the following: 1. Calculate the coordinates of the minimum point,

[0181]

[0182] 2. Initialize the neighborhood tree and insert the neighborhood nodes of the minimum point (mm, nn),

[0183] RBTree→insert({Q (mm-1,nn) , (mm - 1, nn)});

[0184] RBTree→insert({Q (mm+1，nn) , (mm + 1, nn)});

[0185] RBTree→insert({Q (mm,nn-1) , (mm, nn - 1)});

[0186] RBTree→insert({Q (mm，nn+1) , (mm, nn + 1)}).

[0187] 3. Initialize the unwrapping marked graph to 0 and mark the point (mm, nn) as 1: Map_uwp (mm,nn) = 1; 4. Initialize the value of the unwrapped phase diagram to 0 and calculate the phase value of the point (mm, nn): IM_uwp (mm，nn) = P (mm,nn) ; 5. Initialize the set G to an empty set.

[0188] In addition, specifically, the data structure used in step S2A0 is the red - black tree node structure, and the interpretations of its red - black tree node structure and operation functions are as follows:

[0189] ① Define the quality map value Q (m，n) and the coordinates (m, n) as the node on the red - black tree node = (Q (m，n) , (m, n)).

[0190] ② Red - black tree insertion operation: RBTree→insert(node) (insert if the node does not exist, do not operate if it exists).

[0191] ③ Red - black tree deletion operation: RBTree→erase(node) (delete if the node exists, do not operate if it does not exist).

[0192] ④ Red - black tree operation to get the minimum node, compare with the quality map value Q (m，n) as the node value:

[0193] {Q (mi，ni), (ni, ni)} = RBTree->minimum().

[0194] Step S2B0, update the neighborhood tree RBTree, calculate the minimum coordinate (mi, ni), update the set G, and select the neighborhood calculation mode mod.

[0195] Specifically, step S2B0 is as follows:

[0196] 1) Calculate the minimum node of RBTree:

[0197] {Q (mi，ni) , (mi, 1i)} = RBTree->minimum();

[0198] 2) Delete the minimum node of RBTree:

[0199] RBTree->erase({Q (mi，ni) , (mi, ni)});

[0200] 3) Update the reference point set G:

[0201] Calculate the four-neighborhood reference point set G(mi, ni) of the non-boundary point (mi, ni)

[0202] Upper neighborhood reference point, G_up = nozero(Q (mi-1，ni) *Map_wp (mi-1，ni) );

[0203] Lower neighborhood reference point, G_down = nozeo(Q (mi+1，ni) *Map_uwp (mi+1，ni) );

[0204] Left neighborhood reference point, G_left = nozero(Q (mi，ni-1) *Map_uwp (mi，ni-1) );

[0205] Right neighborhood reference point, G_right = nozero(Q (mi，ni+1) *Map_uwp (mi，ni+1) );

[0206] Then the four-neighborhood reference point set, G(mi, i) = {G_up, G_down, G_left, G_right}.

[0207] 4) Select the neighborhood calculation mode:

[0208] ① If the point (mi, ni) is a boundary, that is, mi = 0 or mi = M - 1 or ni = 0 or ni = N - 1, then mod = border.

[0209] ②If mi is not the upper boundary and the reference point in the upper neighborhood is the smallest in the four-neighborhood G(mi, ni), that is, mi ≠ 0 and G_up = min(G(mi, ni)), then mod = up.

[0210] ③If mi is not the lower boundary and the reference point in the lower neighborhood is the smallest in the four-neighborhood G(mi, ni), that is, mi ≠ M - 1 and G_down = min(G(mi, ni)), then mod = down.

[0211] ④If ni is not the left boundary and the reference point in the left neighborhood is the smallest in the four-neighborhood G(mi, ni), that is, ni ≠ 0 and G_left = min(G(mi, ni)), then mod = left.

[0212] ⑤If ni is not the right boundary and the reference point in the right neighborhood is the smallest in the four-neighborhood G(mi, ni), that is, ni ≠ N - 1 and G_right = min(G(mi, ni)), then mod = right.

[0213] Step S2C0: Calculate and unwrap the phase diagram IM_uwp according to the calculation mode classification, and update the neighborhood marking tree RBTree and the unwrapped marking graph Map_uwp.

[0214] Specifically, step S2C0 includes the following steps:

[0215] 1) mod = border

[0216] ①Phase of the unwrapping point (mi, ni): IM_uwp (mi，ni) = 0;

[0217] ②Update the unwrapped marking graph: Map_uwp (mi，ni) = 1;

[0218] 2) mod = up

[0219] ①Phase of the unwrapping point (mi, ni): IM_uwp (mi，ni) = Unwrap(IM_uwp (mi-1，ni) , P (mi，ni) );

[0220] ②Update the unwrapped marking graph: Map_uwp (mi，ni) = 1;

[0221] ③Update the neighborhood tree

[0222] RBTree→insert({Q (mi+1，ni) , (mi + 1, ni)}) if Map_uwp(mi + 1, ni) = 0;

[0223] RBTree -> insert({Q (mi，ni-1) , (mi, ni - 1)}) if Map_uwp(mi, ni - 1) == 0;

[0224] RBTree -> insert({Q (mi，ni+1) , (mi, ni + 1)}) if Map_uwp(mi, ni + 1) == 0.

[0225] 3) mod = down

[0226] ① Unwrap the phase of the unwrapping point (mi, ni): IM_uwp (mi，ni) = Unwrap(IM_uwp (mi+1，ni) , P (mi，ni) );

[0227] ② Update the unwrapped marker map: Map_uwp (mi，ni) = 1;

[0228] ③ Update the neighborhood tree:

[0229] RBTree -> insert({Q (mi-1，ni ), (mi - 1, ni))) if Map_uwp(mi - 1, ni) == 0;

[0230] RBTree -> insert({Q (mi，ni-1) , (mi, ni - 1)}) if Map_uwp(mi, ni - 1) == 0;

[0231] RBTree -> insert({Q (mi，ni+1) , (mi, ni + 1)}) if Map_uwp(mi, ni + 1) == 0.

[0232] 4) mod = left

[0233] ① Unwrap the phase of the unwrapping point (mi, ni): IM_uwp (mi，ni) = Unwrap(IM_uwp (mi，ni-1) , P (mi，ni) );

[0234] ② Update the unwrapped marker map: Map_uwp (mi，ni) = 1;

[0235] ③ Update the neighborhood tree:

[0236] RBTree -> insert({Q (mi-1，ni) , (mi - 1, ni)}) if Map_uwp(mi - 1, ni) == 0;

[0237] RBTree→insert({Q (mi+1，ni ),(mi + 1, ni)}) if Map_uwp(mi + 1, ni) == 0;

[0238] RBTree→insert({Q (mi，ni+1) , (mi, ni + 1)}) if Map_uwp(mi, ni + 1) == 0.

[0239] 5) mod = right

[0240] ① Unwrap the phase of the unwrapping point (mi, ni): IM_uwp (mi，ni) = Unwrap(IM_uwp (mi，ni+1) , P (mi，ni) );

[0241] ② Update the unwrapped marker map: Map_uwp (mini) = 1;

[0242] ③ Update the neighborhood tree:

[0243] RBTree→insert({Q (mi-1，ni) , (mi - 1, ni))) if Map_uwp(mi - 1, ni) == 0;

[0244] RBTree→insert({Q (mi+1，ni) , (mi + 1, ni)}) if Map_uwp(mi + 1, ni) == 0;

[0245] RB Tree→insert({Q (mi，ni-1 ), (mi, ni - 1)}) if Map_uwp(mi, ni - 1) == 0.

[0246] In step S2D0, repeat step S2B0 and step S2C0 until the neighborhood tree RBTree is empty, then all points are unwrapped.

[0247] Specifically, the schematic diagram when all points mentioned in step S2D0 are unwrapped can be referred to Figure 7 , and the schematic diagram before the unwrapping is not completed can be referred to Figure 6 .

[0248] The implementation principle of this embodiment is as follows:

[0249] First, analyze and obtain the most suitable phase guidance map, that is, the phase derivative variance quality map, through the phase derivative variance algorithm based on the MRI image.

[0250] Then, on the basis of obtaining the most suitable phase guidance map, apply the global unwrapping algorithm or the fusion algorithm of the global unwrapping algorithm and the red-black tree structure to start unwrapping from the low phase quality number and extend it to the high-quality region, effectively ensuring the unwrapping of all regions and obtaining a relatively accurate phase map.

[0251] Based on the same inventive concept, an embodiment of the present invention provides an optimized system for phase unwrapping guided by a quality map, including a memory and a processor. The memory stores a program that can be run on the processor to implement any of the Figures 1 to 7 methods.

[0252] An embodiment of the present application also discloses a terminal, including a memory, a processor, and a program stored on the memory and executable on the processor. When the program is loaded and executed by the processor, it can implement any of the Figures 1 to 7 methods.

[0253] The embodiments of the present specific implementation manners are all preferred embodiments of the present application, and do not limit the protection scope of the present application accordingly. Therefore, any equivalent changes made according to the structure, shape, and principle of the present application shall be covered within the protection scope of the present application.

Claims

1. An optimized method for phase unwrapping guided by a quality map, characterized in that, Including: Calculating a phase derivative variance quality map based on an MRI image; Guided by the phase derivative variance quality map, applying a phase unwrapping algorithm to start unwrapping from the point with the lowest phase quality value and transition to the high-quality region until all high-quality regions are unwrapped to obtain a phase map; The phase unwrapping algorithm includes: Initialization: Let the unwrapped phase diagram be a matrix IM_uwp of size M*N, the neighborhood marking diagram be a matrix Map_adj of size M*N, the unwrapped marking diagram be a matrix Map_uwp of size M*N, and the set of four-neighborhood reference points be G, which specifically includes the following: Calculate the coordinates of the minimum value point, Initialize the neighborhood marking diagram to 0 and mark the neighborhood of point (mm, nn) as 1: Map_adj (mm-1，nn) = 1; Map_adj (mm+1，nn) = 1; Map_adj (mm，nn-1) = 1; Map_adj (mm，nn+1) = 1; Initialize the value of the unwrapped point map to 0, mark the point (mm, nn) as 1: Map_uwp (mm，nn) = 1; Initialize the value of the phase map after unwrapping to 0, calculate the phase value of the point (mm, nn): IM_uwp (mm，nn) = P (mm，nn) ; Initialize the set G as an empty set; and, Updating the quality map Q, calculating the minimum value coordinates, updating the set G, and selecting a neighborhood calculation mode mod; Calculating the unwrapped phase map IM_uwp according to the calculation mode classification, and updating the neighborhood marking map Map_adj and the marking map Map_uwp; Repeatedly execute updating the quality map Q, calculating the minimum value coordinates, updating the set G, selecting the neighborhood calculation mode mod, and updating the quality map Q, calculating the minimum value coordinates, updating the set G, selecting the neighborhood calculation mode mod until the neighborhood marking map Map_adj is all 0, then all points are unwrapped; Alternatively, the phase unwrapping algorithm includes: Initialization, define data structures: Let the unwrapped phase diagram be a matrix IM_uwp of size M*N, the neighborhood data structure be a red-black tree RBTree, the marked diagram of unwrapped points be a matrix Map_uwp of size M*N, and the set of four-neighborhood reference points be G, specifically as follows: Calculate the coordinates of the minimum value point, Initialize the neighborhood tree, insert the neighborhood nodes of the minimum value point (mm, nn), RBTree→insert({Q (mm-1,nn) ,(mm - 1, nn)}); RBTree→insert({Q (mm+1,nn) ,(mm + 1,nn)}); RBTree→insert({Q (mm,nn-1) ,(mm,nn - 1)}); RBTree→insert({Q (mm,nn+1) ,(mm,nn + 1)}); Initialize the unwrapped marker graph to 0 and mark the point (mm, nn) as 1: Map_uwp (mm,nn) = 1; Initialize the value of the unwrapped phase graph to 0 and calculate the phase value of the point (mm, nn): IM_uwp (mm,nn) = P (mm,nn) ; Initialize the set G as an empty set; and, Updating the neighborhood tree RBTree, calculating the minimum value coordinates (mi, ni), updating the set G, and selecting a neighborhood calculation mode mod; Calculating the unwrapped phase map IM_uwp according to the calculation mode classification, and updating the neighborhood marking tree RBTree and the unwrapped marking map Map_uwp; Repeatedly execute updating the neighborhood tree RBTree, calculating the minimum value coordinates (mi, ni), updating the set G, selecting the neighborhood calculation mode mod, and calculating the unwrapped phase map IM_uwp according to the calculation mode classification, updating the neighborhood marking tree RBTree and the unwrapped marking map Map_uwp until the neighborhood tree RBTree is empty and all points are unwrapped, where M and N are the height and width of the MRI image respectively, m and n are integer coordinates used to determine the position of pixel points in the image, and mi and ni are variables representing the coordinates of specific pixel points during the calculation process.

2. The optimized method for phase unwrapping guided by a quality map according to claim 1, characterized in that, Calculating the phase derivative variance quality map based on the MRI image includes: An initial phase map is calculated based on the MRI image. Herein, the MRI image is defined as IMG, with a height of M and a width of N. IMG (m,n) represents the value at the m-th row and n-th column of the matrix IMG, where both m and n are integers; the neighborhood set of the point with coordinates (m, n) and its adjacent coordinates above, below, left, and right is set as follows: C{(m,n)}×{(m,n),(m + 1,n),(m - 1,n),(m,n - 1),(m,n + 1)}; Calculating the X-direction phase gradient map; Calculating the y-direction phase gradient map; Calculating the neighborhood mean value of the (m,n) point of the X-direction phase gradient map; Calculating the neighborhood mean value of the (m,n) point of the y-direction phase gradient map; Calculating the phase quality map Q with a height of M and a width of N.

3. The optimized method for phase unwrapping guided by a quality map according to claim 2, characterized in that, The initial phase map calculated based on the MRI image is as follows: P (m，n) = arctan(IMG (m，n) where m ∈ [0, M), n ∈ [0, N), and P (m,n) is the initial phase map.

4. The optimized method for phase unwrapping guided by a quality map according to claim 2, characterized in that, The specific method for calculating the neighborhood mean value of the (m,n) point of the X-direction phase gradient map is as follows: where dx represents the coordinate increment in the x direction when calculating the phase gradient map, and is the average value of the dx values related to the coordinates (m, n).

5. The optimized method for phase unwrapping guided by a quality map according to claim 2, characterized in that, The specific method for calculating the neighborhood mean value of the (m,n) point of the y-direction phase gradient map is as follows: where dy represents the coordinate increment in the y direction when calculating the phase gradient map, is to take the average of the dy values related to the coordinates (m, n).

6. The optimized method for phase unwrapping guided by a quality map according to claim 2, characterized in that,Calculating the phase quality map Q includes: Initialize to the maximum value: Q (m，n) = 100 m ∈ [0, M), n ∈ [0, N); where dx represents the coordinate increment in the x direction when calculating the phase gradient map, is the average value of the dx quantity related to the coordinates (m, n), and dy represents the coordinate increment in the y direction when calculating the phase gradient map, is the average value of the dy quantity related to the coordinates (m, n).

7. An optimization method for phase unwrapping guided by a quality map according to claim 1, wherein, Updating the neighborhood tree RBTree, calculating the minimum value coordinates (mi, ni), and updating the set G, and the selection of the neighborhood calculation mode mod includes: Calculate the minimum node of the RBTree: {Q (mi，m) , (mi, ni)} = RBTree->minimum(); Delete the minimum node of RBTree: RBTree->erase({Q (mi，ni) , (mi, ni)}); Updating the reference point set G: calculating the four-neighborhood reference point set G(mi,ni) of the non-boundary point (mi,ni), and the upper neighborhood reference point G_up = nonzero(Q (mi-1,ni) *Map_uwp (mi-1,ni) ); Lower neighborhood reference point G down=nozero(Q (mi+ni) *Map_uwp (mi+1,ni) );Left neighborhood reference point G_left = nonzero(Q (mi，ni-1) *Map_uwp (mi，ni-1) ); Right neighborhood reference point G_right = nonzero(Q (mi，ni+1) *Map)uwp (mi，ni+1) ); Then the four-neighborhood reference point set G(mi, ni) = {G_up, G_down, G_left, G_right}; Select the neighborhood calculation mode: If the point (mi, ni) is on the boundary, that is, mi = 0 or mi = M - 1 or ni = 0 or ni = N - 1, then mod = border; if mi is not the upper boundary and the upper neighborhood reference point is the smallest in the four-neighborhood G(mi, ni), that is, mi ≠ 0 and G_up = min(G(mi, ni)), then mod = up; if mi is not the lower boundary and the lower neighborhood reference point is the smallest in the four-neighborhood G(mi, ni), that is, mi ≠ M - 1 and G_down = min(G(mi, ni)), then mod = down; if ni is not the left boundary and the left neighborhood reference point is the smallest in the four-neighborhood G(mi, ni), that is, ni ≠ 0 and G_left = min(G(mi, ni)), then mod = left; if ni is not the right boundary and the right neighborhood reference point is the smallest in the four-neighborhood G(mi, ni), that is, ni ≠ N - 1 and G_right = min(G(mi, ni)), then mod = right.

8. An optimization system for phase unwrapping guided by a quality map, wherein, It includes a memory, a processor, and a program stored on the memory and executable on the processor. When the program is loaded and executed by the processor, it can implement an optimization method for phase unwrapping guided by a quality map as described in any one of claims 1 to 7.

Images