A velocity un-wrapping method for four-dimensional blood flow magnetic resonance imaging
By combining image post-processing with constraints of three-dimensional space and time, the velocity winding problem in four-dimensional blood flow magnetic resonance imaging was solved, achieving high signal-to-noise ratio and fast velocity dewinding effect, which is suitable for imaging at different scales such as cardiovascular and cerebrovascular diseases.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- FUDAN UNIVERSITY
- Filing Date
- 2023-02-28
- Publication Date
- 2026-07-10
AI Technical Summary
Current four-dimensional flow magnetic resonance imaging suffers from velocity winding problems, which lead to incorrect flow velocity values. Existing methods often result in low signal-to-noise ratios or excessively long scan times, and existing phase unwinding methods are not suitable for four-dimensional flow magnetic resonance imaging.
A velocity-based dewinding method based on image post-processing is adopted, which combines the three-dimensional spatial continuity and temporal velocity continuity of the blood vessel region of interest. Dewinding is performed through data preprocessing and velocity dewinding steps, including processing based on temporal continuity, with temporal constraints as the main constraint and spatial constraints as the auxiliary constraint, as well as temporal-spatial correlation constraints.
Without increasing the scan time, the signal-to-noise ratio of the velocity distribution map is improved, and the velocity unwinding of four-dimensional blood flow magnetic resonance imaging at different scales, such as cardiovascular and cerebrovascular diseases, can be performed quickly and effectively.
Smart Images
Figure CN115963440B_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of magnetic resonance blood flow imaging technology, specifically relating to a velocity unwinding method for four-dimensional blood flow magnetic resonance imaging. Background Technology
[0002] Four-dimensional flow magnetic resonance imaging (4D MRI) is a dynamic imaging technique that can comprehensively visualize blood vessels in three dimensions and perform hemodynamic analysis. Typically, during imaging, when the actual flow velocity in a certain direction within the imaging range exceeds the pre-set velocity encoding value (VENC) for that direction, phase wrapping occurs in the original image, leading to velocity wrapping in the final flow velocity quantization distribution map. The velocity value at the wrapping point will display an incorrect result. Currently, there are two main solutions to this problem: one is to select an appropriate VENC value before 4D MRI, and the other is to directly de-wrap the velocity based on the obtained velocity distribution map. In the first method, one approach is to set a sufficiently high VENC value based on empirical values of blood flow velocity in the scanned area for subsequent scans; this usually results in a low signal-to-noise ratio in the velocity distribution map. Another approach is to estimate an appropriate VENC value by performing additional 2D pre-scans of the target area; this usually further increases the scan time. In the second category of methods, similar to the phase decoupling required in techniques such as field distribution calculation in magnetic resonance imaging, existing phase decoupling methods used in these techniques are not well applied to four-dimensional blood flow magnetic resonance imaging. Most of the published velocity decoupling methods for blood flow magnetic resonance imaging are used in two-dimensional dynamic blood flow magnetic resonance imaging and mainly employ velocity continuity constraints in the time dimension. There are also phase decoupling methods based on the Laplacian operator specifically for four-dimensional blood flow magnetic resonance imaging, which usually require a long computation time and have only been validated on large arteries such as cardiovascular vessels. Summary of the Invention
[0003] The purpose of this invention is to address the velocity entanglement problem that occurs in traditional four-dimensional blood flow magnetic resonance velocity distribution maps acquired using VENC empirical values, and to provide a velocity dewinding method based on image post-processing.
[0004] The velocity dewinding method for four-dimensional blood flow magnetic resonance imaging based on image post-processing provided by this invention combines the continuity of three-dimensional space and the velocity continuity of time dimension within the region of interest in the four-dimensional blood flow magnetic resonance velocity distribution map for interactive constraints, so as to achieve robust velocity dewinding that can be applied to cardiovascular and cerebrovascular blood flow imaging.
[0005] The velocity decoupling method for four-dimensional blood flow magnetic resonance imaging based on image post-processing provided by this invention mainly consists of two parts: data preprocessing and velocity decoupling. The overall implementation framework is as follows:Figure 1 As shown. In the first part of the data preprocessing, the dynamic three-dimensional spatial velocity distribution map obtained from the four-dimensional blood flow magnetic resonance scan is automatically segmented to obtain the vascular mask map of the region of interest. In the second part of the velocity decoupling process, it is divided into three steps: the first step is velocity decoupling based on temporal continuity constraints; the second step is velocity decoupling based primarily on temporal constraints and secondarily on spatial constraints; and the third step is velocity value correction based on temporal-spatial correlation constraints. Finally, the four-dimensional velocity distribution map after decoupling can be obtained based on the data processing results of the first two parts. Specific explanations are as follows:
[0006] Part 1: Data Preprocessing.
[0007] Four-dimensional flow magnetic resonance imaging (MRI) scans were used to acquire a modulus map (Mag) and dynamic three-dimensional spatial velocity distribution maps (VXt, VYt, and VZt) corresponding to three orthogonal directions. Then, based on these four sets of data, a velocity-weighted modulus map projected as a time-dimensional maximum value (MIP) was calculated. A relatively fast and efficient threshold-based segmentation method was then applied to this MIP map to obtain a 0-1 distributed three-dimensional vascular mask map (VessMask), where pixels with a value of 1 correspond to the region of interest (ROI). Finally, the four-dimensional velocity distribution maps (VXt, VYt, VZt) and the three-dimensional vascular mask map (VessMask) were used as outputs of the data preprocessing section for velocity dewinding in the second part.
[0008] Part Two: Speed Unwinding.
[0009] The velocity unwinding in this section involves processing three different orthogonal directions sequentially. Therefore, the dynamic three-dimensional velocity distribution map (VXt / VYt / VZt) corresponding to each orthogonal direction is processed according to the following steps.
[0010] (1) Velocity unwinding based on time-dimensional continuity constraints:
[0011] For each pixel location with a value of 1 in the 3D vessel mask image VessMask, first determine whether the pixel is a non-isolated point, i.e., whether the sum of VessMask values within a 3x3x3 region centered on the pixel is greater than 1. If this condition is met, it is a non-isolated point. Then, perform the following operations on the pixel location: First, extract the time-dimensional velocity value data of the same location in the 3D velocity distribution map corresponding to the current pixel location to form a velocity-time vector VT. Calculate the difference of this vector in the time dimension and the difference in sign. Then, based on these two results, find the time frame location where the absolute value of the time dimension vector difference is greater than VENC, denoted as indT1, and the time frame location where the sign difference of the time dimension vector is 2, denoted as indT2. Take the intersection of indT1 and indT2, denoted as indWrapT. Then, based on the length of indWrapT, process it according to the following three cases:
[0012] (1.1) If the length of indWrapT is 2, then subtract the sign of the velocity value and the product of 2 VENC from the velocity value corresponding to the time frame between two time frame positions in indWrapT, that is, unwind according to v-sign(v)*2VENC.
[0013] (1.2) If the length of indWrapT is greater than 2, then unwrap the velocity values in indWrapT whose signs are consistent with the velocity signs that account for less than half of the total velocity.
[0014] (1.3) For indWrapT that does not meet the above two conditions, no further unwinding is performed for the time being; if the union of indT1 and indT2 is empty, the position of the pixel is recorded as Pos0; after completing the above judgment and operation, the position of all pixels that have undergone this unwinding process is recorded as Pos1, and the corresponding time frame position after unwinding is Pos1_T.
[0015] (2) Speed unwinding with time-dimensional constraints as the primary constraint and spatial-dimensional constraints as the secondary constraint:
[0016] First, sort the time frames recorded as Pos0 in step (1) according to the absolute value of the speed value. At the same time, sort the time frames recorded as Pos1 and Pos1_T according to the number of times they are unwound. Take the first three time frames in the union of the two as the time frames with the highest probability of speed winding, maxPindT, and take the last three time frames in the union of the two as the time frames with the lowest probability of speed winding, minPindT.
[0017] Then, for each pixel position with a value of 1 in the three-dimensional blood vessel mask image VessMask, if the point is not isolated and is not contained in both Pos0 and Pos1, the following operations are performed: First, extract the time dimension velocity value data of the same position in the three-dimensional velocity distribution map corresponding to the current pixel position to form a velocity-time vector VT, and obtain indWrapT according to the method in step (1); then process each pixel position with a length of 1 indWrapT as follows:
[0018] (2.1) Record the time frame position indTsign1 where the sign of the velocity value in the velocity time vector VT is the same as the sign of the velocity value in time frame indWrapT, and record the time frame position indTsign2 where the sign of the velocity value in VT is the same as the sign of the velocity value in time frame indWrapT+1.
[0019] (2.2) If the length of indTsign1 accounts for less than 40% of the total length of all time frames, find all time frames where the absolute value of the velocity value at indTsign1 is greater than 0.2*VENC, and denot them as indT_mod. Determine whether these time frames indT_mod are included in minPindT. If they are not included, unwrap the velocity value of time frame indWrapT and record the updated velocity time series as VT_mod.
[0020] (2.3) If the length of indTsign2 accounts for less than 40% of the total length of all time frames, find all time frames where the absolute value of the velocity value at indTsign2 is greater than 0.2*VENC, and denot them as indT_mod. Determine whether these time frames indT_mod are included in minPindT. If they are not included, unwrap the velocity value at time frame indWrapT+1 and record the updated velocity time series as VT_mod.
[0021] (2.4) For cases where indT_mod contains minPindT, perform the following operations for each time frame in indT_mod: ① Unwind the velocity value corresponding to the time frame and record the velocity time sequence calculated through unwinding as VT_mod2; ② Summate the absolute value of the difference between the velocity value of the current time frame and the velocity value of the previous time frame in VT_mod with the absolute value of the difference between the velocity value of the current time frame and the velocity value of the next time frame to obtain dist1_T, and sum the absolute values of the difference between the velocity value of the current time frame and the velocity value of the previous time frame in VT_mod2 with the absolute values of the difference between the velocity value of the current time frame and the velocity value of the next time frame. ③ The absolute value of the difference between the velocity value at the pixel in the current time frame of VT_mod and the average velocity value in the adjacent 3x3x3 region in the three-dimensional space of the pixel is recorded as dist1_S, and the absolute value of the difference between the velocity value at the pixel in the current time frame of VT_mod2 and the average velocity value in the adjacent 3x3x3 region in the three-dimensional space of the pixel is recorded as dist2_S; ④ If the ratio of dist1_S to dist2_S is greater than 0.7 or the ratio of dist1_S to dist2_S is greater than 0.5 and dist1_T is greater than dist2_T, then VT_mod2 is used as the updated time series VT_mod at the pixel.
[0022] (2.5) Further determine whether the velocity value sequence of the current pixel should be corrected; if the ratio of the difference in length between indTsign1 and indTsign2 to the total number of time frames is less than 0.5, then record other pixels in the 3x3x3 range of the current pixel that are included in Pos1 or not included in the union of indT1 and indT2. Then, average the velocity time sequence of the pixels that meet these conditions in the direction of the pixel. The correlation between the velocity time vector VT of the current pixel and the average velocity time vector is called coef1, and the correlation between the velocity time vector VT_mod of the current pixel that has just been corrected and the velocity time vector corresponding to the average is called coef2. If coef1>0.5 and coef1>coef2, then record the position of these pixels as Pos2_check, and leave it for further judgment in the third step.
[0023] (3) Velocity value correction based on temporal-spatial correlation constraints:
[0024] First, check the uncertain pixels Pos2_check in step (2); here, each pixel in Pos2_check is processed as follows:
[0025] (3.1) Extract the velocity and time dimension information of the current pixel after the first step of unwinding to form a velocity and time vector VT01, and extract the velocity and time dimension information of the current pixel after the second step of unwinding to form a velocity and time vector VT02.
[0026] (3.2) If the pixel at the current position is included in step1_modPnts among other pixels in the 3x3x3 range, then the time dimension velocity value sequence of these pixels after correction in step (2) is averaged in the direction of the pixel. The correlation between the velocity time vector VT01 of the current pixel and the velocity time vector corresponding to the mean is coef1, and the correlation between the velocity time vector VT02 of the current pixel and the velocity time vector corresponding to the mean is coef2. If coef1>0.5 and coef1>coef2, then the velocity time vector of the pixel at the given position is corrected to VT01.
[0027] Then, for each pixel location with a value of 1 in the 3D blood vessel mask image VessMask, if the point is not isolated and is not contained in both Pos0 and Pos1, the following operations are performed respectively:
[0028] (3.4) Extract the velocity-time dimension information of the current pixel after the previous step to form a velocity-time vector VT1;
[0029] (3.5) Based on VT1, find the time frame position indT1 where the absolute value of the time dimension vector difference is greater than VENC, and the time frame position indT2 where the sign difference of the time dimension vector is 2. Take the union of indT1 and indT2 and denote it as indWrapT.
[0030] (3.6) Determine if indWrapT is not empty. Then compare the velocity value of the pixel at the current position in the neighborhood of the three-dimensional space for each time frame. That is, for each time frame, if the difference between the average velocity value of the other pixels in the 3x3x3 range of the current pixel and the velocity value of the current pixel is greater than 2*VENC, then unwrap the velocity value of the pixel. The velocity time vector after this processing is denoted as VT2.
[0031] (3.7) For the current pixel, the velocity value sequence of the time dimension of other pixels in the 3x3x3 range is averaged in the direction of the pixel. The correlation between the velocity time vector VT1 of the current pixel and the velocity time vector corresponding to the mean is coef1. The correlation between the velocity time vector VT2 of the current pixel after unwinding and the velocity time vector corresponding to the mean is coef2. If coef1>0.5 and coef1>coef2, the velocity time vector of the pixel at the given position is corrected to VT1. Otherwise, it is updated to VT2.
[0032] The velocity distribution map of all pixel locations with a value of 1 in the 3D blood vessel mask image VessMask after the above processing is used as the result after velocity value correction.
[0033] The final four-dimensional velocity distribution map is obtained by unwinding the velocity distribution maps in all three directions through the above three steps, and is used for subsequent visualization and quantitative processing.
[0034] Compared with the prior art, the present invention has the following features and advantages:
[0035] (1) This invention does not require additional scanning with different VENC settings. It only requires setting the VENC according to the empirical value of blood vessel flow velocity in different locations. It ensures a high signal-to-noise ratio of velocity distribution map without increasing the scanning time.
[0036] (2) This invention fully utilizes the time dimension, spatial dimension and spatiotemporal correlation constraints of the data obtained from four-dimensional blood flow magnetic resonance imaging in different velocity directions, and can effectively dewind the velocity distribution map in each direction obtained from four-dimensional blood flow magnetic resonance imaging.
[0037] (3) The present invention can perform efficient speed unwinding of four-dimensional blood flow magnetic resonance imaging of different scales such as cardiovascular and cerebrovascular systems with relatively fast computation processing time. Attached Figure Description
[0038] Figure 1 This is a flowchart illustrating the process framework of the present invention.
[0039] Figure 2 Comparison diagrams of velocity before and after unwinding in the head-to-toe direction are shown. Among them, (a) is the velocity distribution diagram in the head-to-toe direction before unwinding, (b) is the velocity distribution diagram after unwinding using the proposed method, (c) is the waveform diagram of velocity change over time at the internal carotid artery before unwinding, and (d) is the waveform diagram of velocity change over time after unwinding using the proposed method.
[0040] Figure 3The diagrams show the velocity distribution in the left and right directions before and after unwinding. (a) shows the velocity distribution in the left and right directions before unwinding, (b) shows the velocity distribution after unwinding using the proposed method, (c) shows the velocity distribution over time at the right middle cerebral artery before unwinding, and (d) shows the velocity distribution over time after unwinding using the proposed method. Detailed Implementation
[0041] The invention will be further described below with reference to specific examples and accompanying drawings.
[0042] Example 1
[0043] This invention utilizes a Philips 3.0T MRI scanner (Ingenia, Philips Healthcare, Best, The Netherlands) with a 32-channel head coil to perform four-dimensional blood flow MRI scans of the subject's cerebral arterial circle. Specific MRI scan parameters are as follows: FOV = 211mm × 211mm × 30mm, spatial resolution of 1mm isovoxel, flip angle of 10°, repetition time (TR) of 5.6ms, echo time (TE) of 3.2ms, velocity encoding (VENC) in all three directions of 80cm / s, and 18 frames per second.
[0044] Based on the results of four-dimensional blood flow magnetic resonance scanning using the above parameters, this invention first performs data preprocessing. Specifically, the process involves calculating the velocity-weighted modulus map (Mag) and the dynamic three-dimensional spatial velocity distribution maps (VXt, VYt, VZt) corresponding to the three orthogonal directions, respectively. The maximum value projection (MIP) of the modulus map along the time dimension is then calculated. Next, a threshold-based segmentation method is used to obtain a 0-1 distributed three-dimensional vascular mask map (VessMask) from this MIP map; here, the threshold is set to 0.1. Therefore, the four-dimensional velocity distribution maps (VXt, VYt, VZt) and the preprocessed three-dimensional vascular mask map (VessMask) are used together for the next step of velocity dewinding.
[0045] In the subsequent velocity unwinding process, the dynamic three-dimensional velocity distribution map (VXt / VYt / VZt) corresponding to each orthogonal direction is processed according to the following steps:
[0046] (1) Velocity unwinding based on time-dimensional continuity constraints:
[0047] For each pixel location with a value of 1 in the 3D vessel mask image VessMask, first determine if the pixel is a non-isolated point, i.e., whether the sum of VessMask values within a 3x3x3 region centered on the pixel is greater than 1. If this condition is met, it is a non-isolated point. Then, the following operations can be performed on the pixel location: First, extract the time-dimensional velocity value data of the same location in the 3D velocity distribution map corresponding to the current pixel location to form a velocity-time vector VT. Calculate the difference of this vector in the time dimension and the difference in sign. Then, based on these two results, find the time frame location where the absolute value of the time dimension vector difference is greater than 80 (VENC value), denoted as indT1, and the time frame location where the sign difference of the time dimension vector is 2, denoted as indT2. Take the intersection of indT1 and indT2, denoted as indWrapT. Then, based on the length of indWrapT, process it in three ways as follows:
[0048] 1) If the length of indWrapT is 2, then subtract the sign of the velocity value and the product of 2 times VENC from the velocity value corresponding to the time frame between two time frame positions in indWrapT, that is, unwind according to v-sign(v)*2*80.
[0049] 2) If the length of indWrapT is greater than 2, then unwrap the velocity values in indWrapT whose signs are consistent with the velocity signs that account for less than half of the total velocity.
[0050] 3) For indWrapT that does not meet either of the above two conditions, no further unwinding is performed. If the union of indT1 and indT2 is empty, the position of that pixel is recorded as Pos0. After completing the above judgment and operation, all pixel positions that have undergone this unwinding process are recorded as Pos1, and the corresponding time frame position after unwinding is Pos1_T.
[0051] (2) Speed unwinding with time-dimensional constraints as the primary constraint and spatial-dimensional constraints as the secondary constraint:
[0052] First, sort the time frames recorded as Pos0 in step one according to the absolute value of the velocity value. At the same time, sort the time frames recorded as Pos1 and Pos1_T according to the number of unwinding operations. Take the first three time frames in the union of the two as the time frames with the highest probability of velocity winding, maxPindT, and take the last three time frames in the union of the two as the time frames with the lowest probability of velocity winding, minPindT.
[0053] Then, for each pixel position with a value of 1 in the 3D vessel mask image VessMask, if the point is not isolated and is not contained in both Pos0 and Pos1, the following operation is performed. First, extract the time dimension velocity value data of the same position in the 3D velocity distribution map corresponding to the current pixel position to form a velocity-time vector VT, and obtain indWrapT according to the method in step (1). Then, process each pixel position with a length of 1 indWrapT as follows:
[0054] 1) Record the time frame position indTsign1 where the sign of the velocity value in the velocity time vector VT is the same as the sign of the velocity value in time frame indWrapT, and record the time frame position indTsign2 where the sign of the velocity value in VT is the same as the sign of the velocity value in time frame indWrapT+1.
[0055] 2) If the length of indTsign1 accounts for less than 40% of the total length of all time frames, find all time frames where the absolute value of the velocity value at indTsign1 is greater than 0.2*80, and denot them as indT_mod. Determine whether these time frames indT_mod are contained in minPindT. If they are not contained in any of them, unwrap the velocity value of time frame indWrapT and record the updated velocity time series as VT_mod.
[0056] 3) If the length of indTsign2 accounts for less than 40% of the total length of all time frames, find all time frames where the absolute value of the velocity value at indTsign2 is greater than 0.2*80, and denot them as indT_mod. Determine whether these time frames indT_mod are contained in minPindT. If they are not contained in any of them, unwrap the velocity value at time frame indWrapT+1 and record the updated velocity time series as VT_mod.
[0057] 4) For cases where indT_mod contains elements within minPindT, perform the following operations for each time frame in indT_mod: ① Unwind the velocity value corresponding to that time frame and record the velocity time sequence calculated through unwinding as VT_mod2; ② Summate the absolute value of the difference between the velocity value of the current time frame and the velocity value of the previous time frame in VT_mod with the absolute value of the difference between the velocity value of the current time frame and the velocity value of the next time frame to obtain dist1_T. Summate the absolute value of the difference between the velocity value of the current time frame and the velocity value of the previous time frame in VT_mod2 with the absolute value of the difference between the velocity value of the current time frame and the velocity value of the next time frame to obtain dist1_T. ③ The absolute value of the difference between the velocity value at the pixel in the current time frame of VT_mod and the average velocity value in the adjacent 3x3x3 region in the three-dimensional space of the pixel is denoted as dist1_S, and the absolute value of the difference between the velocity value at the pixel in the current time frame of VT_mod2 and the average velocity value in the adjacent 3x3x3 region in the three-dimensional space of the pixel is denoted as dist2_S; ④ If the ratio of dist1_S to dist2_S is greater than 0.7 or the ratio of dist1_S to dist2_S is greater than 0.5 and dist1_T is greater than dist2_T, then VT_mod2 is used as the updated time series VT_mod at the pixel;
[0058] 5) Further determine whether the velocity value sequence of the current pixel should be corrected. If the ratio of the difference in length between indTsign1 and indTsign2 to the total number of time frames is less than 0.5, then record other pixels within the 3x3x3 range of the current pixel that are included in Pos1 or not included in the union of indT1 and indT2. Then, average the velocity time sequences of the pixels that meet these conditions in the pixel direction. The correlation between the velocity time vector VT of the current pixel and the average velocity time vector is recorded as coef1, and the correlation between the velocity time vector VT_mod of the current pixel after correction and the velocity time vector corresponding to the average is recorded as coef2. If coef1>0.5 and coef1>coef2, then record the position of these pixels as Pos2_check, for further judgment in the subsequent third step.
[0059] (3) Velocity value correction based on temporal-spatial correlation constraints:
[0060] First, check the uncertain pixels Pos2_check from step (2). Each pixel in Pos2_check is processed as follows:
[0061] 1) Extract the velocity and time dimension information of the current pixel after the first step of unwinding to form a velocity and time vector VT01, and extract the velocity and time dimension information of the current pixel after the second step of unwinding to form a velocity and time vector VT02.
[0062] 2) If the pixel at the current position is included in step1_modPnts among other pixels within a 3x3x3 range, then the time dimension velocity value sequence of these pixels after the second correction is averaged along the pixel direction. The correlation between the velocity time vector VT01 of the current pixel and the velocity time vector corresponding to the average value is defined as coef1, and the correlation between the velocity time vector VT02 of the current pixel and the velocity time vector corresponding to the average value is defined as coef2. If coef1>0.5 and coef1>coef2, then the velocity time vector of the pixel at the given position is corrected to VT01.
[0063] Then, for each pixel location with a value of 1 in the 3D blood vessel mask image VessMask, if the point is not isolated and is not contained in both Pos0 and Pos1, the following operations are performed respectively:
[0064] 1) Extract the velocity-time dimension information of the current pixel after the previous step to form a velocity-time vector VT1;
[0065] 2) Based on VT1, find the time frame position indT1 where the absolute value of the time dimension vector difference is greater than 80, and the time frame position indT2 where the sign difference of the time dimension vector is 2. Take the union of indT1 and indT2 and denote it as indWrapT.
[0066] 3) Check indWrapT. If indWrapT is not empty, compare the velocity value of the pixel at the current position in the neighborhood of the three-dimensional space for each time frame. That is, for each time frame, if the difference between the average velocity value of other pixels in the 3x3x3 range of the current pixel and the velocity value of the current pixel is greater than 2*80, then unwrap the velocity value of the pixel. The velocity time vector after this processing is denoted as VT2.
[0067] 4) For the current pixel, the time dimension velocity value sequence of other pixels within a 3x3x3 range is averaged along the pixel direction. The correlation between the current pixel's velocity time vector VT1 and the velocity time vector corresponding to the average value is defined as coef1. The correlation between the current pixel's velocity time vector VT2 (which has just been unwound) and the velocity time vector corresponding to the average value is defined as coef2. If coef1 > 0.5 and coef1 > coef2, then the velocity time vector of the pixel at the given position is corrected to VT1; otherwise, it is updated to VT2.
[0068] The velocity distribution map of all pixel locations with a value of 1 in the 3D blood vessel mask image VessMask after the above processing is used as the result after velocity value correction.
[0069] The final four-dimensional velocity distribution map is obtained by unwinding the velocity distribution maps in all three directions through the above three steps, and is used for subsequent visualization and quantitative processing.
[0070] Below are some examples of results.
[0071] Figure 2 This paper presents a comparison of the velocity distribution map along the head-to-foot direction of one layer obtained from a four-dimensional blood flow magnetic resonance imaging (MRI) scan before and after unwinding. The comparison of the velocity distribution maps and the velocity-time waveforms before and after unwinding in this image demonstrates that the proposed velocity unwinding method can effectively recover pixels with velocity winding from the original scan.
[0072] Figure 3 This paper presents a comparison of the velocity distribution map along the left-right direction of one layer obtained from a four-dimensional blood flow magnetic resonance imaging (MRI) scan before and after unwinding. The image shows that the original acquisition results exhibited large-scale, even discontinuous, velocity winding at the middle cerebral arteries on both sides. The proposed velocity unwinding method effectively restores these velocity-wound pixels.
Claims
1. A velocity unwinding method for four-dimensional blood flow magnetic resonance imaging, characterized in that, By combining the three-dimensional spatial continuity and temporal velocity continuity characteristics within the blood vessel region of interest in the four-dimensional blood flow magnetic resonance velocity distribution map, interactive constraints are applied to achieve robust velocity decoupling that is applicable to both cardiovascular and cerebrovascular blood flow imaging; specifically, it consists of two parts: data preprocessing and velocity decoupling. Data preprocessing involves automatically segmenting the dynamic three-dimensional spatial velocity distribution map obtained from a four-dimensional blood flow magnetic resonance scan to obtain a vascular mask map of the region of interest. Speed unwinding involves three steps: (1) Velocity unwinding based on time-dimensional continuity constraints; (2) For speed unwinding, the time dimension is the main constraint and the space dimension is the auxiliary constraint; (3) The velocity value is corrected based on the temporal and spatial correlation constraints; finally, the four-dimensional velocity distribution map after unwinding is obtained based on the data processing results of the first two parts. The data preprocessing specifically includes: First, a modulus map (Mag) and dynamic three-dimensional spatial velocity distribution maps (VXt, VYt, and VZt) corresponding to three orthogonal directions are obtained through four-dimensional blood flow magnetic resonance scanning. Then, based on these four sets of data, a velocity-weighted modulus map is calculated, which is a projection MIP map of the maximum value along the time dimension. Next, a threshold-based segmentation method is used on this MIP map to obtain a 0-1 distributed three-dimensional vascular mask map (VessMask), where the position of a pixel with a value of 1 corresponds to the vascular region of interest (ROI). Finally, the four-dimensional velocity distribution maps (VXt, VYt, VZt) and the three-dimensional vascular mask map (VessMask) are used as the output of the data preprocessing part for velocity dewinding in the second part. In the velocity unwinding step (1), for each pixel with a value of 1 in the three-dimensional blood vessel mask image VessMask, it is first determined whether the pixel is not an isolated point, that is, whether the sum of the values of VessMask in the 3x3x3 area centered on the pixel is greater than 1. If the condition is met, it is a non-isolated point. Then, the following operations are performed on the pixel: First, extract the time dimension velocity value data at the same position in the three-dimensional velocity distribution map corresponding to the current pixel to form a velocity time vector VT. Calculate the difference of the vector in the time dimension and the difference of the positive and negative signs. Then, based on these two results, find the time frame position where the absolute value of the time dimension vector difference is greater than VENC and record it as indT1, and the time frame position where the sign difference of the time dimension vector is 2 and record it as indT2. Take the intersection of indT1 and indT2 and record it as indWrapT. VENC is the set velocity encoding value. In the step (2) of velocity unwinding, firstly, the time frames recorded as Pos0 in step (1) are sorted according to the absolute value of the velocity value, and the time frames recorded as Pos1 and Pos1_T are sorted according to the number of unwinding times. The first three time frames in the union of the two are taken as the time frame maxPindT with the highest probability of velocity winding, and the last three time frames in the union of the two are taken as the time frame minPindT with the lowest probability of velocity winding. Then, for each pixel with a value of 1 in the three-dimensional blood vessel mask image VessMask, if the point is not an isolated point and is not included in Pos0 and Pos1 at the same time, the following operation is performed: firstly, the time dimension velocity value data of the same position in the three-dimensional velocity distribution map corresponding to the current pixel point is extracted to form a velocity time vector VT, and indWrapT is obtained according to the method in step (1). In step (3) of speed unwinding, firstly, check the uncertain pixel points Pos2_check from step (2); process each pixel point in Pos2_check as follows: (3.1) Extract the velocity and time dimension information of the current pixel after unwinding in step (1) to form a velocity and time vector VT01, and extract the velocity and time dimension information of the current pixel after unwinding in step (2) to form a velocity and time vector VT02; (3.2) For other pixels within the 3x3x3 range of the current pixel, if Pos1 is included, the time dimension velocity value sequence of these pixels after correction in step (2) is averaged in the direction of the pixel. The correlation between the velocity time vector VT01 of the current pixel and the velocity time vector corresponding to the mean is coef1, and the correlation between the velocity time vector VT02 of the current pixel and the velocity time vector corresponding to the mean is coef2. If coef1>0.5 and coef1>coef2, the velocity time vector of the pixel at the given position is corrected to VT01.
2. The velocity unwinding method for four-dimensional blood flow magnetic resonance imaging according to claim 1, characterized in that, In the step (1) of speed unwinding, the following three cases are handled according to the length of indWrapT: (1.1) If the length of indWrapT is 2, then subtract the sign of the velocity value and the product of 2 VENC from the velocity value corresponding to the time frame between two time frame positions in indWrapT, that is, unwind according to v-sign(v)*2VENC; (1.2) If the length of indWrapT is greater than 2, then unwrap the velocity values in indWrapT whose signs are consistent with the velocity signs that account for less than half of the total velocity. (1.3) For indWrapT that does not meet the above two conditions, no further unwinding is performed for the time being; if the union of indT1 and indT2 is empty, the position of the pixel is recorded as Pos0; after completing the above judgment and operation, the position of all pixels that have undergone this unwinding process is recorded as Pos1, and the corresponding time frame position after unwinding is Pos1_T.
3. The velocity unwinding method for four-dimensional blood flow magnetic resonance imaging according to claim 2, characterized in that, In step (2) of the speed unwinding, each pixel with a length of 1 indWrapT is processed as follows: (2.1) Record the time frame position indTsign1 where the sign of the velocity value in the velocity time vector VT is the same as the sign of the velocity value in time frame indWrapT, and record the time frame position indTsign2 where the sign of the velocity value in VT is the same as the sign of the velocity value in time frame indWrapT+1. (2.2) If the length of indTsign1 accounts for less than 40% of the total length of all time frames, find all time frames where the absolute value of the velocity value at indTsign1 is greater than 0.2*VENC, and record them as indT_mod. Determine whether these time frames indT_mod are included in minPindT. If they are not included, unwrap the velocity value of time frame indWrapT and record the updated velocity time series as VT_mod. (2.3) If the length of indTsign2 accounts for less than 40% of the total length of all time frames, find all time frames where the absolute value of the velocity value at indTsign2 is greater than 0.2*VENC, and record them as indT_mod. Determine whether these time frames indT_mod are included in minPindT. If they are not included, unwrap the velocity value at time frame indWrapT+1 and record the updated velocity time series as VT_mod. (2.4) For cases where indT_mod contains minPindT, the following operations are performed on each time frame in indT_mod: ① Unwind the velocity value corresponding to the time frame and record the velocity time sequence calculated by unwinding as VT_mod2; ② Sum the absolute values of the difference between the velocity value of the current time frame and the velocity value of the previous time frame in VT_mod and the absolute values of the difference between the velocity value of the current time frame and the velocity value of the next time frame to obtain dist1_T; sum the absolute values of the difference between the velocity value of the current time frame and the velocity value of the previous time frame in VT_mod2 and the absolute values of the difference between the velocity value of the current time frame and the velocity value of the next time frame to obtain dist2_T; ③ The absolute value of the difference between the velocity value at the pixel in the current time frame of VT_mod and the average velocity value of the adjacent 3x3x3 region in the three-dimensional space of the pixel is denoted as dist1_S. The absolute value of the difference between the velocity value at the pixel in the current time frame of VT_mod2 and the average velocity value of the adjacent 3x3x3 region in the three-dimensional space of the pixel is denoted as dist2_S. If the ratio of dist1_S to dist2_S is greater than 0.7 or the ratio of dist1_S to dist2_S is greater than 0.5 and dist1_T is greater than dist2_T, then VT_mod2 is used as the updated time series VT_mod at the pixel. (2.5) Further determine whether the velocity value sequence of the current pixel should be corrected; if the ratio of the difference between the lengths of indTsign1 and indTsign2 to the total number of time frames is less than 0.5, then record the other pixels in the 3x3x3 range of the current pixel if they are included in Pos1 or not included in the union of indT1 and indT2, then average the velocity time sequence of the pixels that meet these conditions in the direction of the pixel, and respectively set the correlation between the velocity time vector VT of the current pixel and the mean velocity time vector as coef1, and set the correlation between the velocity time vector VT_mod of the current pixel that has just been corrected and the velocity time vector corresponding to the mean as coef2. If coef1>0.5 and coef1>coef2, then record the position of these pixels as Pos2_check, and leave it for further judgment in subsequent step (3).
4. The velocity unwinding method for four-dimensional blood flow magnetic resonance imaging according to claim 3, characterized in that, In step (3) of speed unwinding, for each pixel with a value of 1 in the three-dimensional blood vessel mask image VessMask, if the point is a non-isolated point and is not contained in Pos0 and Pos1 at the same time, the following operations are performed respectively: (3.4) Extract the velocity-time dimension information of the current pixel after the previous step to form a velocity-time vector VT1; (3.5) Based on VT1, find the time frame position indT1 where the absolute value of the time dimension vector difference is greater than VENC, and the time frame position indT2 where the sign difference of the time dimension vector is 2. Take the union of indT1 and indT2 and denote it as indWrapT. (3.6) Judge indWrapT. If indWrapT is not empty, compare the velocity value of the pixel at the current position in the neighborhood of the three-dimensional space for each time frame. That is, for each time frame, if the difference between the average velocity value of other pixels in the 3x3x3 range of the pixel at the current position and the velocity value of the current pixel is greater than 2*VENC, then unwrap the velocity value of the pixel. The velocity time vector after this processing is denoted as VT2. (3.7) For the current pixel, the velocity value sequence of the time dimension of other pixels in the 3x3x3 range is averaged in the direction of the pixel. The correlation between the velocity time vector VT1 of the current pixel and the velocity time vector corresponding to the mean is coef1. The correlation between the velocity time vector VT2 of the current pixel after unwinding and the velocity time vector corresponding to the mean is coef2. If coef1>0.5 and coef1>coef2, the velocity time vector of the pixel at the given position is corrected to VT1. Otherwise, it is updated to VT2. The velocity distribution map of all pixels with a value of 1 in the 3D blood vessel mask image VessMask after the above processing is used as the result after velocity value correction.