Fourier Velocity Encoded MRI: Acceleration and Velocity Map Estimation

Abstract

Fourier velocity encoding (FVE) is an alternative to phase contrast imaging (PC). FVE provides considerably higher SNR than PC, due to its higher dimensionality and larger voxel sizes. Furthermore, FVE is robust to partial voluming, as it resolves the velocity distribution within each voxel. FVE data are usually acquired with low spatial resolution, due to scan-time restrictions associated with its higher dimensionality. FVE is capable of providing the velocity distribution associated with a large voxel, but does not directly provides a velocity map. Knowing the velocity distribution on a voxel is important for accurate diagnosis of stenosis in vessels on the scale of spatial resolution. Velocity maps, however, are useful for visualizing the actual blood flow through a vessel and can be used in different studies and diagnosis. In this context, this chapter deals with two aspects of the FVE MRI technique: acceleration and estimation of velocity map. First, are introduced six different acceleration techniques that can be applied to FVE acquisition. Methods such as variable-density sampling and compressive sampling. Then, is proposed a novel method to estimate velocity maps with high spatial resolution from low-resolution FVE data. Finally, it can be concluded that FVE datasets can be acquired in time scale comparable to PC, it contains more velocity information, since it resolves a velocity distribution within a voxel, and also provides an accurate estimation of the velocity map.

1. Introduction

Cardiovascular diseases are among the main causes of death in both men and women in the United States. Some of these diseases are caused or can be diagnosed by abnormal blood flow in a particular part of the cardiovascular system. For example, atherosclerosis consists of the narrowing of a blood vessel due to the gradual accumulation of lipids, inflammatory cells and connective tissue in the vessel wall [1]. This narrowing alters the local blood flow and may cause flow jets and/or turbulent flow. In these flow jets occur peaks of velocity that are significantly higher than those exhibited at a normal flow. Thus, knowledge of blood flow patterns in the human body is an important component in the research and diagnosis of certain cardiovascular diseases. Currently, two distinct approaches to the study and quantification of blood flow in the human body are available to researchers and clinicians: in-vivodirect measurements of the velocity field using velocity-encoded magnetic resonance imaging (MRI) or Doppler ultrasound.

Doppler ultrasound is the gold standard for quantifying blood flow patterns in the clinical environment. The equipment is relatively small, cheap and portable, and is capable of producing measurements in real time with excellent temporal resolution. On the other hand, evaluation by ultrasound is inadequate when there is fat, air, bone, or surgical scar in the acoustic path. Moreover the equipment is strongly user-dependent, since flow measurements are inaccurate when the ultrasound beam cannot be properly aligned with the axis of flow [2, 3].

MRI is capable of three-dimensional visualization of all aspects of a cardiac examination, such as the anatomy of the heart, features in the blood vessels, and also the quantification of velocity in any given vessel. Compared to ultrasound, magnetic resonance imaging does not have the same operator dependence, being able to accurately quantify the correct direction of flow, and does not have the same acoustic limitations related to bones, fat, air or surgical scars.

The current gold standard for MRI flow quantification is phase contrast (PC) [4]. In this technique, a bipolar gradient is aligned to the flow axis to obtain a velocity measurement (approximately the mean [5]) for each voxel of the image. Despite its unrestricted use, phase contrast has some limitations. Phase contrast technique suffers from partial-volume effects when a wide distribution of velocities is contained within a single voxel [6]. This is particularly problematic when flow is turbulent and/or complex (e.g., flow jets due to stenosis) or at the interface between blood and vessel wall (viscous sublayer). This issue is typically addressed by increasing the spatial resolution, which dramatically affects the signal-to-noise ratio (SNR) and increases the scan time. Therefore, PC may be inadequate for estimating the peak velocity of stenotic flow jets and for assessing wall shear rate.

Fourier velocity encoded (FVE) MRI [7] is a magnetic resonance velocity quantification technique which is as an alternative to phase contrast imaging, since real-time FVE is the MRI equivalent to spectral-Doppler ultrasound [8]. In this technique, the acquired measurements have a considerably higher signal-to-noise ratio than those acquired with phase contrast, due to its high-dimensional data set and also to its larger voxels. In addition, different from PC data, FVE does not suffer from partial volume effects, since for each voxel a velocity distribution is measured. So this technique can accurately diagnose vessels stenosis on low spatial resolution. The data set measured with this technique is usually obtained with very low spatial resolution. This is due to restrictions associated with its high dimensionality, which can lead to long acquisitions time. Thus, FVE is not a popular technique in the clinical environment that requires exams to be performed as fast as possible. On the other hand, it has been shown that the FVE acquisition can be accelerated. For example, FVE acquisition using rapid spiral sampling in $k$-space is a fast and reliable alternative to accurately measure velocity peaks in blood flow jets or to obtain hemodynamic parameters [9].

In this context, this chapter deals with two aspects of the FVE MRI technique: acceleration and estimation of velocity map. First, are introduced six different important acceleration techniques that can be applied to FVE acquisition and are related to the use of variable-density sampling, which may be used along spatial k-space and velocity k-space, partial Fourier acquisition along velocity k-space, temporal acceleration methods such as UNFOLD and k-t BLAST, parallel imaging methods and compressive sampling.

Finally, since FVE does not provide the actual velocity map associated with the flow, is proposed a novel method to velocity maps estimation with high spatial resolution from low-resolution FVE data. The proposed method is based on the mathematical model of the FVE distribution, $s\left(x,y,v\right)$, and involves solving a PDE-constrained optimization related to the Navier-Stokes equation.

2. Magnetic resonance flow imaging

MRI is a modality uniquely capable of imaging all aspects of cardiovascular disease, and is a potential “one-stop shop” for cardiovascular health assessment. MRI can generate cross-sectional images in any plane (including oblique planes), and can also measure blood flow. The image acquisition is based on using strong magnetic fields and non-ionizing radiation in the radio frequency range, which are harmless to the patient. MR is used to image hydrogen nuclei, because of its abundance in the human body. Spinning charged particles (or “spins”), such as hydrogen nuclei, act like a tiny bar magnet, presenting a very small magnetic field, emanating from the south pole to the north pole. In this section we introduce the mathematical formalism of MR imaging and flow imaging.

2.1. Mathematical formalism

The acquired MR signal $s\left(t\right)$at a particular time instant corresponds to a sample of the Fourier transform $M\left(k_x,k_y\right)$of the excited magnetization $m\left(x,y\right)$:

$$M\left(k_x,k_y\right)={\int }_x{\int }_{y}m\left(x,y\right)e^{-j2\pi \left(k_{x}x+k_{y}y\right)}\mathit{dx}\mathit{dy}.$$E1

The Fourier coordinates $k_x$and $k_y$vary with time, according to the zeroth moment of the readout gradients $G_x$and $G_y$:

$$k_x\left(t\right)=\frac{\gamma }{2\pi }{\int }_0^t{G}_x\left(\tau \right)\mathit{d\tau }$$E2

$$k_y\left(t\right)=\frac{\gamma }{2\pi }{\int }_0^t{G}_y\left(\tau \right)\mathit{d\tau }.$$E3

This formalism can be generalized for any combination of $G_x$, $G_y$, and $G_z$gradients:

$$M\left({{\displaystyle \overrightarrow{k}}}_r\right)={\int }_{{\displaystyle \overrightarrow{r}}}m\left(\overrightarrow{r}\right)\cdot e^{-j2\pi {{\displaystyle \overrightarrow{k}}}_r\cdot {\displaystyle \overrightarrow{r}}}d{\displaystyle \overrightarrow{r}}$$E4

$${{\displaystyle \overrightarrow{k}}}_r\left(t\right)=\frac{\gamma }{2\pi }{\int }_0^t{{\displaystyle \overrightarrow{G}}}_r\left(\tau \right)\mathit{d\tau },$$E5

where ${{\displaystyle \overrightarrow{G}}}_r$is the oblique gradient resulting from the combination of the $G_x$, $G_y$and $G_z$gradients, and $\overrightarrow{r}$is its corresponding axis along which the linear variation in magnetic field intensity is realized.

Given a spatial position function ${\displaystyle \overrightarrow{r}}$(t) and a magnetic field gradient ${{\displaystyle \overrightarrow{G}}}_r\left(t\right)$, the magnetization phase is:

$$\varphi \left({\displaystyle \overrightarrow{r}},,,t\right)=\gamma {\int }_0^t{{\displaystyle \overrightarrow{G}}}_r\left(\tau \right)\cdot \overrightarrow{r}\left(\tau \right)\mathit{d\tau },$$E6

For static spins, $\overrightarrow{r}\left(t\right)$is constant ($\overrightarrow{r}$), and this becomes:

$$\varphi =\gamma \overrightarrow{r}\cdot {\int }_0^t{{\displaystyle \overrightarrow{G}}}_r\left(\tau \right)\mathit{d\tau }$$E7

$$=2\pi {{\displaystyle \overrightarrow{k}}}_r\cdot \overrightarrow{r},$$E8

as in the exponential in Eq. (4).

2.2. Principles of MR flow imaging

The basic principles of quantitative flow measurement using magnetic resonance were first proposed by Singer [10] and Hahn [11] in the late 1950s. However, clinical applications of MR flow quantification were not reported until the early 1980s [12, 13, 14, 15]. Current MR flow imaging methods are based on the fact that spins moving at a constant velocity accrue a phase proportional to the velocity times the first moment of the gradient waveform along the direction in which they are moving.

For spins moving along the $\overrightarrow{r}$-axis with a constant velocity $\overrightarrow{v}$, and initial position ${{\displaystyle \overrightarrow{r}}}_0$, we can write $\overrightarrow{r}\left(t\right)={{\displaystyle \overrightarrow{r}}}_0+{\displaystyle \overrightarrow{v}}t$. Rewriting Eq. (6), for $t=t_0$:

$$\varphi =\gamma {\int }_0^{t_0}{{\displaystyle \overrightarrow{G}}}_r\left(t\right)\cdot \left({{\displaystyle \overrightarrow{r}}}_0+{\displaystyle \overrightarrow{v}}t\right)\mathit{dt}$$E9

$$=\gamma {{\displaystyle \overrightarrow{r}}}_0\cdot {\int }_0^{t_0}{{\displaystyle \overrightarrow{G}}}_r\left(t\right)\mathit{dt}+\gamma {\displaystyle \overrightarrow{v}}\cdot {\int }_0^{t_0}{{\displaystyle \overrightarrow{G}}}_r\left(t\right)t\mathit{dt}$$E10

$$=\gamma {{\displaystyle \overrightarrow{r}}}_0\cdot {{\displaystyle \overrightarrow{M}}}_0+\gamma {\displaystyle \overrightarrow{v}}\cdot {{\displaystyle \overrightarrow{M}}}_1,$$E11

where ${{\displaystyle \overrightarrow{M}}}_0$and ${{\displaystyle \overrightarrow{M}}}_1$are the zeroth and first moments of the $\overrightarrow{r}$-gradient waveform at the time of signal acquisitions (“echo time”, or “time to echo” (TE)), respectively. Thus, if a gradient with null zeroth moment is used (e.g., a bipolar gradient, aligned with $\overrightarrow{v}$), the phase accrued for a constant velocity spin is $\varphi =\gamma {\displaystyle \overrightarrow{v}}\cdot {{\displaystyle \overrightarrow{M}}}_1$.

Therefore, if a bipolar gradient waveform is played between the excitation and the readout, the phase measured in a pixel of the acquired image is directly proportional to the velocity of the spins contained within its corresponding voxel. However, factors other than flow (such as inhomogeneities of the magnetic field) may cause additional phase shifts that would cause erroneous interpretation of the local velocity [16].

2.2.1. Phase contrast

The phase contrast method addresses the problem mentioned above by using two gradient-echo data acquisitions in which the first moment of the bipolar gradient waveform is varied between measurements [4]. So from Eq. (11) it is possible to obtain time-dependent velocity measures in all three spatial directions. Then for a fixed time and direction, e.g. velocity in $z$-axis, the through-plane velocity in each voxel is measured as:

$$v_z\left(x,y\right)=\frac{{\varphi }_a\left(x,y\right)-{\varphi }_b\left(x,y\right)}{\gamma \left(M_1^a-M_1^b\right)},$$E12

where ${\varphi }_a\left(x,y\right)$and ${\varphi }_b\left(x,y\right)$are the phase images acquired in each acquisition, and $M_1^a$and $M_1^b$are the first moment of the bipolar gradients used in each acquisition.

2.2.2. Fourier velocity encoding

While phase contrast provides a single velocity measurement associated with each voxel, Fourier velocity encoding [7] provides a velocity histogram for each spatial location, which is a measurement of the velocity distribution within each voxel.

FVE involves phase-encoding along a velocity dimension. Instead of only two acquisitions, as in phase contrast, multiple acquisitions are performed, and a bipolar gradient with a different amplitude (and first moment) is used in each acquisition. Eq. (11) can be rewritten as:

$$\varphi \left(\overrightarrow{r},\overrightarrow{v},t\right)=2\pi \left({{\displaystyle \overrightarrow{k}}}_r\cdot {\displaystyle \overrightarrow{r}}+{{\displaystyle \overrightarrow{k}}}_v\cdot {\displaystyle \overrightarrow{v}}\right),$$E13

where ${{\displaystyle \overrightarrow{k}}}_v$is the velocity frequency variable associated with $\overrightarrow{v}$, and is proportional to the first moment of ${{\displaystyle \overrightarrow{G}}}_r\left(t\right)$:

$${{\displaystyle \overrightarrow{k}}}_v=\frac{\gamma }{2\pi }{{\displaystyle \overrightarrow{M}}}_1.$$E14

Each voxel of the two-dimensional image is associated with a distribution of velocities. This three-dimensional function $m\left(x,y,v\right)$is associated with a three-dimensional Fourier space $M\left(k_x,k_y,k_v\right)$. Thus, an extra dimension is added to k-space, and multiple acquisitions are required to cover the entire $k_x$-$k_y$-$k_v$space ( Figure 1 ). In order to move along $k_v$, a bipolar gradient with the appropriate amplitude (and first moment) is played before the $k_x$-$k_y$readout gradients, in each acquisition. Placing the bipolar gradient along the $z$-axis will encode through-plane velocities. Placing the bipolar gradient along $x$or $y$will encode in-plane velocities. Oblique flow can be encoded using a combination of bipolar gradients along the $x$, $y$and $z$axes.

Figure 1.

Spiral FVE k-space sampling scheme. The dataset corresponding to each temporal frame is a stack-of-spirals in$k_x‐k_y‐k_v$space. Each spiral acquisition corresponds to a different$k_v$encode level.

Each acquisition along $k_v$is called a velocity encode. The number of required velocity encodes depends on the desired velocity resolution and velocity field-of-view (the maximum range of velocities measured without aliasing). For example, to obtain a 25 cm/s resolution over a 600 cm/s field-of-view, 24 velocity encodes are needed. The spatial-velocity distribution, $m\left(x,y,v\right)$, is obtained by inverse Fourier transforming the acquired data, $M\left(k_x,k_y,k_v\right)$. If cine imaging [17] is used, measurements are also time resolved, resulting in a four-dimensional dataset: $m\left(x,y,v,t\right)$.

2.3. FVE signal model

2DFT phase contrast provides two 2-dimensional functions, $m\left(x,y\right)$and $v_z\left(x,y\right)$, the magnitude and velocity maps, respectively. For simplicity we are assuming that the through-plane velocity map is in the $z$direction. If these maps are measured with sufficiently high spatial resolution, and flow is laminar, one can assume that each voxel contains only one velocity, and therefore the spatial-velocity distribution associated with the object is approximately:

$$s\left(x,y,v\right)=m\left(x,y\right)\times \delta \left(v-v_z\left(x,y\right)\right),$$E15

where $\delta \left(v\right)$is the Dirac delta function.

In 2DFT FVE, k-space data is truncated to a rectangular cuboid in $k_x$-$k_y$-$k_v$space. The associated object domain spatial-velocity blurring can be modeled as a convolution of the true object distribution, $s\left(x,y,v\right)$, with $\mathit{sinc}\left(x/\mathrm{\Delta }x\right)$, $\mathit{sinc}\left(y/\mathrm{\Delta }y\right)$, and $\mathit{sinc}\left(v/\mathrm{\Delta }v\right)$, where $\mathrm{\Delta }x$and $\mathrm{\Delta }y$are the spatial resolutions along the $x$and $y$axes, respectively, and $\mathrm{\Delta }v$is the velocity resolution, as follows:

$$\widehat{s}\left(x,y,v\right)=\left[m\left(x,y\right)\times \delta (v-v_z(x,y\left)\right)\right]\ast sinc\left(\frac{x}{\mathrm{\Delta }x}\right)\ast sinc\left(\frac{y}{\mathrm{\Delta }y}\right)\ast sinc\left(\frac{v}{\mathrm{\Delta }v}\right),$$E16

where $\widehat{s}\left(x,y,v\right)$is the measured object distribution and $\ast$denotes convolution. This is equivalent to:

$$\widehat{s}\left(x,y,v\right)=\left[m\left(x,y\right)\times sinc\left(\frac{v-v_z\left(x,y\right)}{\mathrm{\Delta }v}\right)\right]\ast \left[sinc\left(\frac{x}{\mathrm{\Delta }x}\right)\times sinc\left(\frac{y}{\mathrm{\Delta }y}\right)\right].$$E17

On the other hand, spiral FVE acquisitions follows a stack-of-spirals pattern in $k_x$-$k_y$-$k_v$space ( Figure 1 ), then k-space data is truncated to a cylinder, i.e., a circle along $k_x$-$k_y$(with diameter $1/\mathrm{\Delta }r$), and a rect function along $k_v$(with width $1/\mathrm{\Delta }v$), where $\mathrm{\Delta }r$and $\mathrm{\Delta }v$are the prescribed spatial and velocity resolutions, respectively. Using the same approach we used for 2DFT FVE, the associated object domain spatial-velocity blurring in spiral FVE can be modeled as a convolution of the true object distribution, $s\left(x,y,v\right)$, with $\mathrm{\text{jinc}}\left(\sqrt{x^2+y^2}/\mathrm{\Delta }r\right)$and $\mathit{sinc}\left(v/\mathrm{\Delta }v\right)$, resulting in:

$$\begin{array}{l}\widehat{s}\left(x,y,v\right)=\left[m\left(x,y\right)\times \delta (v-v_z(x,y\left)\right)\right]\ast sinc\left(\frac{v}{\mathrm{\Delta }v}\right)\ast \text{jinc}\left(\frac{\sqrt{x^2+y^2}}{\mathrm{\Delta }r}\right)\\ & =\left[m\left(x,y\right)\times sinc\left(\frac{v-v_z\left(x,y\right)}{\mathrm{\Delta }v}\right)\right]\ast \text{jinc}\left(\frac{\sqrt{x^2+y^2}}{\mathrm{\Delta }r}\right),\hfill \end{array}$$E18

where $\text{jinc}\left(z\right)=J_1\left(\mathit{\pi z}\right)/\left(2z\right)$and $J_1\left(z\right)$is the Bessel function of the first kind and first order. These approaches for deriving FVE data from high-resolution velocity maps will be used for the map estimation purposes.

3. Acceleration of FVE

FVE datasets are multidimensional, which makes this method particularly suitable for accelerated acquisition. Variable-density sampling may be used along spatial k-space, and also along velocity k-space. Partial Fourier acquisition along velocity k-space can be used to reduce scan time by nearly 50%. Temporal acceleration methods such as UNFOLD and k-t BLAST have been demonstrated with FVE. Parallel imaging methods have also been shown to work well with FVE. Also FVE is optimally suited for acquisition acceleration using compressed sensing. This section introduces each of these acceleration methods.

3.1. Variable-density sampling of spatial k-space

Magnetic resonance imaging can be accelerated using variable-density sampling of k-space. This is typically implemented by using a sampling pattern that satisfies the Nyquist criterion at the low spatial frequencies, and undersamples the high spatial frequencies. In other words, the effective field-of-view (FOV) is varied from the desired FOV at the center of k-space to a reduced FOV at the periphery [18]. The general hypothesis is that artifacts from undersampling the periphery of k-space will be negligible, because the energy of high frequency components is typically much lower than that of low frequency components. Variable-density spirals can increase spatiotemporal resolution and improve accuracy in flow quantitation [19]. The spatial aliasing resulting from variable-density spiral sampling is incoherent, and, in the regions-of-interest (e.g., cardiac chambers, valves, great vessels), it typically originates from static or slow moving material located at the periphery of the spatial FOV (e.g., chest wall). FVE resolves the distribution of velocities within the voxel, thus moderate low-velocity aliasing artifacts generally do not affect one’s ability to calculate diagnostically important parameters—such as peak velocity and acceleration—from the time-velocity distribution.

The use of variable-density spirals for acceleration of slice-selective FVE with spiral readouts is illustrated in Figure 2 . A single-shot uniform-density spiral readout was replaced with a multi-shot variable-density spiral acquisition. The use of multi-shot acquisitions provides the possibility of multi-dimensional temporal acceleration, and allows reduction of readout duration and TR, which reduce off-resonance artifacts and temporal aliasing, respectively. The use of a shorter TR also allows improving the temporal resolution. The data in Figure 2a was obtained using a single-interleave 8 ms readout uniform-density spiral design [20, 21]. The variable-density design used three 4 ms spiral interleaves, and provided higher spatial resolution and reduced off-resonance artifacts, and thus better spatial localization of flow ( Figure 2b ) [9]. Some aliasing artifacts were observed in spatial domain (see asterisk), but these were not observed in the time-velocity distributions. A fully sampled reference is shown in Figure 2c , for comparison.

Figure 2.

Effect of variable-density sampling of spatial k-space on image quality and spatial localization of flow: (a) uniform-density design; (b) variable-density design; (c) ground truth reference. Top row: spatial images from the first cardiac phase; center row: time-velocity distributions measured at the aortic valve; bottom row: time-velocity distributions measured in the descending aorta. The use of higher spatial resolution and shorter readout duration improves the spatial localization of flow, which is identified by the reduced signal from static material in the time-velocity histograms (see arrows). Some aliasing artifacts were observed in spatial domain (see asterisk), but these were not observed in the time-velocity distributions.

3.2. Variable-density sampling of velocity k-space

Variable-density sampling of velocity k-space was first demonstrated by DiCarlo et al. [22] using real-time FVE. Real-time FVE (also known as MR Doppler or one-shot FVE) [8, 23, 24, 25] utilizes cylindrical excitation to restrict the spatial field-of-view to a one-dimensional beam. An oscillating readout gradient simultaneously encodes spatial position and velocity along the axis of the beam. Variable-density sampling of velocity k-space has also been demonstrated using slice-selective FVE [26]. Variable-density sampling along the velocity dimension may be used to improve the velocity resolution and/or increase the velocity field-of-view. However, conventional non-Cartesian reconstruction methods such as gridding and direct Fourier transform (DrFT) do not adequately deal with the associated undersampling artifacts. Alternatively, reconstruction of variable-density FVE may be performed using variable-width sinc interpolation with dynamic field-of-view centering [26]. Figure 3 illustrates the use of variable-density sampling along velocity k-space for accelerating slice-selective FVE [26]. The reconstruction scheme using variable-width sinc interpolation with dynamic field-of-view centering exhibits negligible aliasing artifacts compared to conventional gridding (see arrows). There is also no noticeable loss of velocity resolution compared with the small velocity FOV ground truth reference. Note the improvement in velocity resolution compared with the large FOV uniform-density result.

Figure 3.

In vivodemonstration of variable-density sampling of velocity k-space. Velocity distributions were measured using slice-selective spiral FVE at the aortic valve plane of a healthy volunteer using: (a) uniform-density sampling, large FOV; (b) uniform-density sampling, small FOV (ground truth); (c) variable-density sampling, reconstructed using conventional gridding; and (d) variable-density sampling, reconstructed using variable-width sinc interpolation with dynamic field-of-view centering. The reconstruction scheme using variable-width sinc interpolation with dynamic field-of-view centering reduces undersampling artifacts (arrows), and shows velocity resolution equivalent to that of the ground truth reference.

3.3. Partial Fourier acquisition of velocity k-space

Partial Fourier acquisition and reconstruction exploits the conjugate symmetry property of the Fourier transform of real-valued signals. The method involves acquiring slightly greater than one half of k-space, and synthesizing the missing data using a combination of conjugate synthesis and background phase correction. A narrow strip of k-space is acquired with symmetric coverage in order to estimate this smoothly-varying background phase. The fastest and most widely used method of partial Fourier reconstruction is homodyne detection [27]. Acquisition time in FVE can be reduced by 30–40% using partial Fourier acceleration along the velocity dimension. This consists in acquiring only slightly more than half of the $k_v$encodings, and synthesizing the missing data using homodyne reconstruction. This has been successfully used in FVE for scan time reduction, without significant loss of velocity resolution. This approach has been demonstrated in studies with healthy volunteers [8, 20, 21] and patients [20, 21, 25, 28], and in phantom experiments [22]. The feasibility of reducing scan time in FVE using partial Fourier acquisition is illustrated in Figure 4 . Up to 42% of the acquired data (along velocity k-space) was discarded and then synthesized using homodyne reconstruction. The results show 71 and 60% improvement in velocity resolution using this approach, when imaging the aortic valves of a healthy volunteer and of a patient with aortic stenosis. Partial Fourier performs well in both healthy volunteer and patient studies, and no significant loss of resolution or artifacts is noticed [20, 21].

Figure 4.

Evaluation of partial k-space reconstruction along the velocity dimension, in aortic valve studies of a healthy volunteer (a–c) and a patient with aortic stenosis (d–f). Homodyne reconstruction performs well in both healthy volunteer (b) and patient (e) studies, improving the velocity resolution by 71 and 60%, respectively. Full k-space distributions with the same number of velocity-encode samples are shown for comparison (a,d), as well as the fully sampled datasets (c,f).

3.4. Temporal acceleration

In dynamic MRI, view sharing [29] is commonly used to increase the number of temporal frames. Artifacts and loss of temporal resolution due to view sharing can be avoided or corrected using temporal acceleration techniques, such as UNFOLD [30, 31] and k-t BLAST [32]. UNFOLD reduces scan time by making efficient use of k-t space, and can be very successful in the context of slice-selective FVE due to the high dimensionality of this imaging method. The use of UNFOLD for acceleration of FVE was first demonstrated by Macgowan and Madore [33], and further investigated by Carvalho and Nayak [9, 20, 21]. Figure 5 illustrates an implementation of the UNFOLD method specially designed for slice-selective FVE with spiral readouts [9, 20, 21]. A view-ordering scheme that reduces overlap in v-f space was designed (v denotes the through-plane velocity dimension, and f denotes temporal frequency). Figure 2a shows the undersampled data in both $v$–$f$and $v$–$t$domains (where t denotes time). The aliasing signal is filtered using a two-dimensional filter ( Figure 5a ). This filter has a bandwidth of 107 Hz for velocities below $\pm$150 cm/s. For higher velocities, the bandwidth varies from 69 to 30 Hz. This results in effective temporal resolutions of 9.3 and 14.5–33.3 ms, respectively. The temporal resolution is lower for higher velocities, but this may prove unnoticeable, as the velocity distribution of high-velocity flow jets within large voxels is typically temporally smooth. For comparison, the temporal resolution with view sharing would be 50 ms for all velocities ( Figure 5d ). The remaining narrow-bandwidth aliasing components at $\pm$20 and $\pm$40 Hz are filtered using a tight zero-phase one-dimensional notch filter along the temporal dimension ( Figure 5b ). The final results show that this temporal acceleration scheme is capable of achieving 6-fold acceleration in multi-interleaf spiral FVE, without noticeable loss of temporal resolution, and without introducing significant artifacts ( Figure 5c ). View-sharing ( Figure 5d ), on the other hand, is equivalent to a moving-average low-pass filter, which reduces the temporal frequency bandwidth (dashed arrows), and causes loss of temporal resolution, perceived as blurring along time (circled).

Figure 5.

Temporal acceleration compared with view sharing in (left) v-f space and (right) v-t space: (a) undersampled data; (b) with two-dimensional filtering; (c) with two-dimensional and notch filtering; and (d) with view sharing. The two-dimensional filter (dashed lines) removes a majority of the aliasing, and the notch filter (dotted line) removes the remaining aliasing signal (solid arrows). This approach removes aliasing components without noticeable loss of temporal resolution. View sharing reduces the temporal frequency bandwidth (dashed arrows) and causes temporal blurring (circles).

3.5. Parallel imaging

Spatial aliasing due to undersampling of slice-selective FVE can be reduced using parallel imaging methods such as SENSE [34] and SPIRiT [35]. Parallel imaging is an acceleration approach that uses data from multiple coils to reduce aliasing artifacts due to undersampling of spatial k-space [34]. Steeden et al. was able to accelerate slice-selective spiral FVE by a factor of four using SENSE [28]. Lyra-Leite et al. used two-dimensional and three-dimensional SPIRiT to accelerate slice-selective spiral FVE by factors of two and four, respectively [36, 37]. In the velocity distributions measured using slice-selective FVE, aliasing due to spatial undersampling typically results in increased signal at $v=0$cm/s, since the majority of the aliasing signal is associated with static material. Figure 6 illustrates the use of two-dimensional SPIRiT to accelerate slice-selective spiral FVE by a factor of two [36]. SPIRiT is able to considerably reduce aliasing artifacts, while not introducing significant artifacts (see error images).

Figure 6.

Time-velocity distributions from select voxels, reconstructed using twofold accelerated two-dimensional SPIRiT (center row), in comparison with the fully sampled reference (top row): (a) right external carotid artery; (b) right internal carotid artery; and (c) left carotid bifurcation.

3.6. Compressive sensing

Compressive sensing (CS) has been used in MRI [38] context for a while in different applications, such as fMRI images [39], PC-MRI velocity maps [40] and also FVE distributions [41, 42]. Basically, is a set of theories and methods that establish the conditions under which a signal can be reconstructed based on a limited number of linear measurements. It also states different procedures for signal reconstruction, provided that these conditions are properly met [43, 44, 45, 46]. For a successful image reconstruction using CS the desired image must satisfy three conditions: (1) must have a sparse representation in a known transform domain, (2) artifacts caused by k-space undersampling must be incoherent in the sparsifying transform domain and (3) must be reconstructed by a nonlinear method that enforces both sparsity of the image representation and consistency of the reconstruction with the acquired samples [38].

FVE data is suitable for CS application, since the information contained in images with different velocity encodes is highly redundant differing only where flow occurs. Therefore, through spatial finite differencing operations FVE dataset have a sparse representation [38, 42].

The original CS reconstruction problem is a NP-hard problem, generally of combinatorial complexity [46, 47, 48], and is not viable except for very low-dimensional cases. Thus, the original problem can be relaxed and a precise reconstruction can be achieved using the following non-linear constrained optimization problem:

$${\mathbf{f}}^{\ast }={\text{argmin}}_{\mathbf{f}}\parallel \mathbf{Tf}{\parallel }_{{\mathrm{\ell }}_p}^p\mathrm{s}.\mathrm{t}.\mathbf{Mf}=\mathbf{s},$$E19

where $0<p\le 1$, $\mathbf{T}$is the sparsifying transform, $\mathbf{M}$is the acquisition process matrix, $\mathbf{f}$is the desired image, $\mathbf{s}$is the acquired signal and

$$\parallel \mathbf{f}{\parallel }_{{\mathrm{\ell }}_p}={\left(\sum _{n=1}^N{\left|f_n\right|}^p\right)}^{1/p}.$$E20

Usually in most CS applications the value of $p$is set to $p=1$, but it has been shown that for ${\mathrm{\ell }}_p$-minimization (with $0<p<1$) requires fewer measurements than ${\mathrm{\ell }}_1$[46]. In order to reconstruct MR data based on ${\mathrm{\ell }}_p$-minimization, one can use the algorithm described by Miosso et al. [45].

Other possible ways to enhance signal reconstruction in CS, both in terms of reducing the number of required measurements and in terms of improving image quality for a fixed number of measurements, include the use of support prior information extracted from structural knowledge, previous frames or previous slices [39, 46], and the use of information extracted using machine learning techniques [49, 50]. Other alternative optimization problems are also desired in the context of noisy measurements, in which case, for example, the equality constraint in Problem 19 is replaced by an inequality such as $\parallel \mathbf{Mf}-\mathbf{s}{\parallel }_{{\mathrm{\ell }}_2}\le \epsilon$, with $\epsilon$being a tolerance to noise [47, 48] — the higher the value of $\epsilon$, the higher the number of measurements required for reconstruction.

In this context, has been shown by Marinelli et al. [51] and Hilbert et al. [42] that CS can also be used as an acceleration technique for FVE datasets and the acquisition can be made in time scale comparable to the gold standard phase contrast. So it is possible to obtain meaningful velocity spectra in small vessels in clinical time while regular phase contrast can provide only mean velocity maps [42].

4. Estimating velocity maps from FVE distributions

In this section will be discussed a methodology to estimate the velocity map based on the FVE velocity distribution. It has been shown in Section 2.3 that FVE velocity distribution signal model $\widehat{s}\left(x,y,v\right)$is related to the actual velocity map $v_z\left(x,y\right)$through the relation

$$\widehat{s}\left(x,y,v\right)=\left[m\left(x,y\right)\times sinc\left(\frac{v-v_z\left(x,y\right)}{\mathrm{\Delta }v}\right)\right]\ast \Psi \left(x,y\right),$$E21

where $\Psi \left(x,y\right)$is a point spread function associated with $k$-space truncation data. This provide a first relation between the FVE measured velocity distribution and the velocity map. On the other hand, blood can be ideally modeled as an incompressible Newtonian fluid. Then, blood flow can be predicted using the Navier-Stokes equation

$$\rho \left(\frac{\partial \mathbf{v}}{\partial t}+\mathbf{v}\cdot \nabla \mathbf{v}\right)=-\nabla p+\mu {\nabla }^2\mathbf{v},$$E22

where $\mathbf{v}=\left(v_x,v_y,v_z\right)$is the velocity vector, $\rho$is the blood density, $\mu$is the whole blood viscosity and ${\nabla }^2$is the Laplacian differential operator. Then, ideally the desired velocity map must satisfy the flow physics model. Therefore, for a fixed instant of time, a velocity map can be estimated from a measured FVE dataset $f\left(x,y,v\right)$, with $K$velocity encodes, through the following PDE-constrained optimization problem

$$\underset{v_z}{min}\sum _{k=1}^K{\int }_{\mathrm{\Omega }}{\left\{f\left(\mathbf{x},v_k\right)-\left[m\left(\mathbf{x}\right)\times sinc\left(\frac{v_k-v_z}{\mathrm{\Delta }v}\right)\right]\ast \Psi \left(\mathbf{x}\right)\right\}}^2\mathit{dA}\mathrm{s}.\mathrm{t}.\rho \mathbf{v}\cdot \nabla \mathbf{v}=-\nabla p+\mu {\nabla }^2\mathbf{v},$$E23

where $\mathbf{x}=\left(x,y\right)$is the position vector and $v_k$is a velocity encode.

In order to solve Eq. (23) the Navier-Stokes equation must be discretized. Since the interest here is in a proof-of-concept velocity map estimation based on only one component of the velocity vector, a bidimensional version of the physics model solver was used. Fluid is assumed incompressible, so the steady 2D Navier-Stokes-continuity dimensionless system of equations [52],

$$\mathbf{v}\cdot \nabla \mathbf{v}=-\nabla p+\frac{1}{\mathit{Re}}{\nabla }^2\mathbf{v}\text{and}\nabla \cdot \mathbf{v}=0,$$E24

was discretized using the Finite Element Method [53], where $\mathit{Re}$is the Reynolds number [52], $\mathbf{v}=v_x\mathbf{i}+v_z\mathbf{j}\in {\mathrm{IR}}^2$is the velocity field and $p$is the pressure. Discretization is made using residues functions based on the governing equations’ weak form Gresho and Sani [53]

$$R_c\left(\mathbf{v}\right)={\int }_{\mathrm{\Omega }}\left(\nabla \cdot \mathbf{v}\right)\mathit{\varphi d}\mathrm{\Omega }$$E25

and

$$R_m\left(\mathbf{v},p\right)={\int }_{\mathrm{\Omega }}\left(\mathbf{v}\cdot \nabla \mathbf{v}\right)\cdot \mathit{\Psi d}\mathrm{\Omega }+{\int }_{\mathrm{\Omega }}\mathrm{s}:\nabla \mathit{\Psi d}\mathrm{\Omega }-{\int }_{\mathrm{\Gamma }}\left(\mathbf{n}\cdot \mathrm{s}\right)\cdot \mathit{\Psi d}\mathrm{\Gamma },$$E26

where $\varphi \in \mathrm{IR},\Psi \in {\mathrm{IR}}^2$are test functions, and $\sigma =-p\mathbf{I}+{\mathit{Re}}^{-1}\left[\nabla \mathbf{v}+\nabla {\mathbf{v}}^T\right]$the Newtonian stress tensor [52].

Discretizatized equations are written as a linear system $\mathbf{Jc}=\mathbf{r}$, where $\mathbf{J}$is a matrix given by the residues’ Jacobian, $\mathbf{r}$is a vector given by the residues and $\mathbf{c}$is the solution vector containing velocity and pressure. Now the minimization problem given by Eq. (23) can be written as

$$\underset{{\mathbf{v}}_z}{min}\sum _{k=1}^K{‖{\mathbf{f}}_k-\left[\mathbf{m}\times sinc\left(\frac{v_k-{\mathbf{v}}_z}{\mathrm{\Delta }v}\right)\right]\ast \Psi ‖}_{{\mathrm{\ell }}_2}^2+\lambda {‖\mathbf{J}\left[{\mathbf{v}}_x;{\mathbf{v}}_z;\mathbf{p}\right]-\mathbf{r}‖}_{{\mathrm{\ell }}_2}^2,$$E27

where $\mathbf{c}=\left[{\mathbf{v}}_x;{\mathbf{v}}_z;\mathbf{p}\right]$is the solution vector written in a stacked form and $\mathbf{m}$is a spin density map with high spatial resolution.

In order to validate the proposed constrained optimization (Eq. (27)) an simple experiment was carried out. To do so, a FVE dataset was simulated from an acquired PC dataset, then the optimization was solved and finally the resultant velocity map was compared with the acquired PC velocity map qualitatively and quantitatively.

First, high-spatial-resolution four-dimensional PC data of a pulsatile carotid flow phantom (Phantoms by Design, Inc., Bothell, WA) were obtained using a 3DFT SPGR pulse sequence. The scan parameters were: $0.5\times 0.5\times 1$mm${}^3$spatial resolution; field-of-view $4.0\times 3.5\times 5.0$cm${}^3$; TR 11.4 ms; flip angle 8.5°; temporal resolution 91.2 ms; VENC 50 cm/s; 40 min per scan; 9 NEX. The data were acquired on a GE Discovery MR750 3T system, with a 32-channel receive-only head coil array (Nova Medical, Inc., Wilmington, MA, USA). The through-slab ($z$) axis was oriented along the S/I direction. The phantom’s pulse cycle was set to 60 bpm. The velocity map for each spatial axis—$u_{\mathrm{pc}},$$v_{\mathrm{pc}}$, and $w_{\mathrm{pc}}$—was reconstructed using data from all channels of the receive coil array. The lumen was segmented by manually outlining the vessel borders from a stack of 2D axial images, obtained from the reconstructed 3D volume.

Then simulated spiral FVE distributions were derived from the acquired phase contrast data using the signal model presented in Eq. (21). Simulated data was generated only for the through-axis velocity component ($v_z$), and for a cardiac phase corresponding to the phantom’s mid-systole. The 9-NEX PC dataset was used in this process, so that the FVE distributions were computed from low-noise velocity maps (as in Carvalho et al. [54]). This is because FVE has considerably higher SNR than PC in general, due to its higher dimensionality and larger voxel size. Finally, two different spiral FVE distributions were obtained for each slice of the volume with $\mathrm{\Delta }r$= 2 mm spatial resolution: one using the proposed method and the other one using the method proposed by Rispoli and Carvalho [55]. The velocity resolution was set to $\mathrm{\Delta }v=10$cm/s, over a 120 cm/s velocity field-of-view.

About the discretization of the Navier-Stokes equations, lumen manually outlined was used to define computational mesh and simulation grid was designed with 1.0 $\times$0.5 ${\mathrm{mm}}^2$element resolution using $Q_2/P_{-1}$elements. Phantom’s blood-mimicking fluid (with Reynolds number $\mathit{\text{Re}}=110$) was assumed to be Newtonian and incompressible. PC-MRI velocity profile was set at the inlet together with no-slip boundary condition.

The optimization problem given by Eq. (27) was then solved using a alternating minimization technique [56]. Left side part was solved using a standard non-linear least squares algorithm and the physics model part of the optimization was solved using Newton’s method [53].

Figure 7 presents the results of the validation experiment using the phase contrast velocity map acquired at the pulsatile carotid flow phantom’s bifurcation. The velocity maps estimated from the simulated low spatial resolution FVE data are very similar (qualitatively) to the reference map. At first glance one can say that the velocity map obtained using the technique proposed by Rispoli and Carvalho [55] ( Figure 7c ) is more similar to the acquired PC-MRI velocity map. However the error images show that the velocity map obtained using the technique proposed in this work ( Figure 7b ) was more accurate than the one obtained with the other method ( Figure 7c ).

Figure 7.

Validation experiment using a pulsatile carotid flow phantom: (a) reference phase contrast velocity map, measured at the phantom’s bifurcation; (b) velocity map estimated from the simulated low-resolution spiral FVE data with$\mathrm{\Delta }r$= 2 mm spatial resolution with the proposed method (and associated error percentages); and (c) velocity map estimated from the simulated low-resolution spiral FVE data with$\mathrm{\Delta }r$= 2 mm spatial resolution with the method proposed by Rispoli and Carvalho [55] (and associated error percentages).

Moreover, a quantitative comparison was performed based on the signal-to-error ratio (SER). The acquired phase contrast velocity field, $v_{\mathrm{pc}}$, was used as the ground-truth “signal”; consequently, the estimation error is the difference between the estimated velocity field, $v_e$, and the ground-truth field, $v_{\mathrm{pc}}$. Thus, the SER is calculated (in decibels) as:

$${\mathrm{SER}}_{{\displaystyle \overrightarrow{\nu }}}=10{log}_{10}\left(\frac{\sum _{i,j}{‖v_{\mathrm{pc}}\left(i,j\right)‖}^2}{\sum _{i,j}{‖v_e\left(i,j\right)-v_{\mathrm{pc}}(i,j)‖}^2}\right),$$E28

Finally, the proposed method measured SER, relative to the PC reference, was 44.63 dB while the technique proposed in Rispoli and Carvalho [55] achieved 28.68 dB. Showing that the proposed optimization given by Eq. (27) is more consistent with the actual velocity map than the previous method proposed.

These good results are important, meaning that FVE may potentially be a substitute of PC imaging, since it contains both a velocity distribution and also velocity map with considerably higher SNR and robustness to partial voluming.

5. Conclusion

In this chapter, was discussed approaches in order to make Fourier Velocity Encoding MRI more suitable for the clinical environment. FVE is a promising MRI technique capable of measuring blood flow in the blood vessels and estimating important biomarkers that are useful for understand and diagnose diseases. It provides a velocity distribution within a voxel instead of a mean velocity map like phase contrast but requires acceleration to be feasible in the clinical setting. So was discussed six different strategies that can reduce drastically the acquisition time. The acceleration techniques discussed are related to the use of variable-density sampling, which may be used along spatial k-space and velocity k-space, partial Fourier acquisition along velocity k-space, temporal acceleration methods such as UNFOLD and k-t BLAST, parallel imaging methods and compressive sensing.

On the other hand, was proposed a novel method for estimating high-resolution velocity maps from low-resolution FVE measurements. This method is based on a PDE-constrained optimization that incorporates the FVE signal model and the Navier-Stokes equation. Results showed that it is possible to obtain highly accurate velocity maps from the FVE distributions. Finally, it can be concluded that FVE datasets can be acquired in time scale comparable to the gold standard phase contrast, it provides more velocity information, since it contains a velocity distribution, and also can provide the actual velocity map as long as a constrained-optimization problem to restore the velocity map is solved.