4D ultrafast ultrasound flow imaging: in vivo quantification of arterial volumetric flow rate in a single heartbeat

We present herein 4D ultrafast ultrasound flow imaging, a novel ultrasound-based volumetric imaging technique for the quantitative mapping of blood flow. Complete volumetric blood flow distribution imaging was achieved through 2D tilted plane-wave insonification, 2D multi-angle cross-beam beamforming, and 3D vector Doppler velocity components estimation by least-squares fitting. 4D ultrafast ultrasound flow imaging was performed in large volumetric fields of view at very high volume rate (>4000 volumes s−1) using a 1024-channel 4D ultrafast ultrasound scanner and a 2D matrix-array transducer. The precision of the technique was evaluated in vitro by using 3D velocity vector maps to estimate volumetric flow rates in a vessel phantom. Volumetric Flow rate errors of less than 5% were found when volumetric flow rates and peak velocities were respectively less than 360 ml min−1 and 100 cm s−1. The average volumetric flow rate error increased to 18.3% when volumetric flow rates and peak velocities were up to 490 ml min−1 and 1.3 m s−1, respectively. The in vivo feasibility of the technique was shown in the carotid arteries of two healthy volunteers. The 3D blood flow velocity distribution was assessed during one cardiac cycle in a full volume and it was used to quantify volumetric flow rates (375  ±  57 ml min−1 and 275  ±  43 ml min−1). Finally, the formation of 3D vortices at the carotid artery bifurcation was imaged at high volume rates.


Introduction
The 3D visualization and quantification of complex blood flow patterns in the human cardiovascular system is a challenging imaging problem that requires the estimation of blood flow velocities along all directions with high spatial and high temporal resolutions within a large field of view. The recent introduction of 4D flow cardiovascular magnetic resonance (CMR) imaging has demonstrated the clinical potential of complete flow analyses that include the quantification of the three components of the velocity flow vectors (Harloff et al 2009, Markl et al 2012. 4D CMR is a powerful yet expensive tool that can image and quantify blood flows at any location in the body with a composite temporal resolution of, typically, 30 ms, by combining hundreds of heart cycles acquired over tens of minutes. Optical coherence tomography is also a promising blood flow imaging technique but remains limited to superficial tissues application (Makita et al 2008, Srinivasan et al 2010, e.g. the human eye vascularization. Conventional ultrasound Doppler imaging is used in daily clinical practice to assess in realtime deep and complex blood flow, with sub-millimeter resolution but it is limited to the estimation of the flow velocity component along the beam axis. Several ultrasound-based vector flow imaging techniques have been developed over the past decades to quantify and visualize both the axial and lateral components of the blood velocity (Bohs et al 2000, Capineri et al 2002, Jensen 2003, Ebbini 2006, Tortoli et al 2006, Udesen and Jensen 2006, Xu and Bashford 2013, Jensen 2014) in a single heartbeat. Real-time implementations of these approaches have been evaluated in vivo and are now available in clinical scanners (Hansen et al 2015, Tanaka et al 2015. Such techniques were also extended to 3D vector flow imaging, but were limited to a small number of imaging planes and the frame-rates achieved remained relatively low (see Dunmire et al 2000 for a review). To achieve higher frame-rates, sequences based on the emission of unfocused waves that insonify large fields of view can be used. For instance, 2D ultrafast ultrasound imaging based on plane-wave emissions has been proposed for high-frame-rate axial velocity (Sandrin et al 1999) or velocity vector (Tanter et al 2002) estimation, in which the increased frame-rate can either be leveraged to map low velocity blood flow with high sensitivity (Macé et al 2011) or to detect high velocities in large arteries (Bercoff et al 2011). Ultrafast ultrasound imaging has been applied to the vector estimation of cardiac (Lovstakken et al 2006) and arterial (Flynn et al 2011, Ekroll et al 2013, Yiu et al 2014, Fadnes et al 2015 blood flow. 2D methods based on diverging-wave emissions, which provide very large field of views, were also developed for the mapping of blood flow in real-time, both for a single (Tanaka et al 2008, Osmanski et al 2012 and for both components of the velocity vector (Pihl et al 2012).
The development of standard 3D ultrasound systems opened new possibilities for vector flow imaging in 3-dimensions by assessing new and reliable measurements, such as 3D velocity vector components and the volumetric flow rate. For instance, volumetric echocardiographic particle image velocimetry (V-Echo-PIV) (Falahatpisheh et al 2014, Falahatpisheh andKheradvar 2016) was shown to assess and quantify in vivo cardiac vector flow in 3D dimensions at up to 30 imaging volumes per second after the injection of microbubbles.
Likewise, the advent of 3D ultrasound systems containing a large number of programmable channels has recently opened new possibilities in ultrafast imaging (Tanter and Fink 2014) and blood flow imaging (Jensen et al 2013). For instance, we recently introduced 3D ultrafast Doppler imaging (Provost et al 2014(Provost et al , 2015, which allowed for the mapping of one component of the blood velocity in large 3D fields of view, e.g. the entire heart (Provost 2014) or the carotid bifurcation (Provost 2015) at thousands of volumes per second. In parallel, the transverse oscillation approach was developed and applied to the estimation of 3D velocity vectors in vitro using a 2D matrix-array probe .
Holbek et al presented an in vivo carotid artery 3D velocity vector study in one plane using continuous data and a 3D transverse oscillations approach (Holbek et al 2015). However, the full 3D mapping of the three components of the velocity vector in a full volumetric field of view at high frame rate remains an important technical challenge.
In this study, we introduce 4D ultrafast flow imaging based on ultrafast cross-beam vector Doppler imaging for the volumetric mapping of blood flow at high volume-rate. The objectives of this study were (1) to establish the experimental methods for volumetric blood flow imaging, (2) to evaluate the accuracy and precision of the technique in vitro in artery phantoms and (3) to show the in vivo feasibility of mapping arterial flows and assessing volumetric flow rates in the human carotid artery.

Multi-angle volumetric ultrafast Doppler acquisition concept
In order to reconstruct the three components of the local velocity vector in a volumetric field of view at thousands of volume per second, 4D ultrafast flow imaging is based on (1) multiple 2D tilted plane-wave insonifications, (2) receive beamforming performed along multiple receive angles (figure 1), (3) spatiotemporal clutter filtering, and (4) voxel-wise least-square velocity vector estimation. ). The emission of the complete set of 2D tilted plane waves was then repeated to obtain a continuous dataset. -angle-delay-and-sum beamforming. Beamforming was performed using a multi-angle-delay-and-sum approach. Specifically, pre-defined receive angles determined an associated voxel-dependent sub-aperture that was then used in a standard delay-and-sum beamforming algorithm. This voxel-wise calculation was repeated for each transmit, for each receive angle, and for each volume.

2D multi
A total number of M receive sub-apertures were defined according to two parameters: the f/D ratio, i.e. the ratio of the imaging depth to the aperture size, and the receive angle . Given the very large number of possible combinations, M was chosen according to computational considerations. The resulting output of the beamforming algorithm was a matrix of size N × M × nx × ny × nz × nt, where nx, ny, and nz are the number of voxels in the x-, y-, and z-dimensions and nt the number of volumes.
2.1.3. Clutter filtering. Spatiotemporal clutter filtering  was applied individually to each M × N data pair. This clutter filtering technique uses the singular value decomposition of the column-wise form obtained from beamformed demodulated data to remove clutter artifacts according to the highly discriminative spatiotemporal coherence of tissue compared to blow flow signals.

Velocity vector and volumetric flow rate estimation.
A phase-shift matrix ( ) Φ IR was calculated using an autocorrelation estimator (Kasai et al 1985), for which the following relationship holds for one ( ) p q x y z t IQ p q x y z t IQ p q x y z t t , , , , , arg , , , , , , , , , , ,

IR PRF
(1) IQ is the 4D beamformed and demodulated data, x and y are the lateral directions, z is the depth, t is the instantaneous time, t PRF is the time between each transmitted pulse and = … p N 1, , and = … q M 1, , . When multiple ( )  →  → k k , i p r q , , pairs are used, equation (1) amounts to a system of equations with more equations than unknowns and it can thus be solved using a least-squares fitting approach via a pseudo-inverse operation. A NM ×3 matrix, K IR , is thus defined as: where indices x, y, and z, indicate the three components of the vectors along the x-, y-, and z-directions. The blood velocity vector, ( ) → v x y z , , , can then be obtained by the Doppler frequency shift equation, as follows: is the pseudo-inverse matrix, c is the speed of sound, f 0 is the transmitted frequency and PRF is the pulse repetition frequency.
The least squares fitting was shown to be robust to noise in 2D implementations (Ekroll et al 2013, Yiu et al 2014. The volumetric flow rate was computed by integrating the flux in a region of interest (ROI), i.e. the rate of flow of a property per unit area. Specifically, after segmenting a ROI using B-mode volumes, the dot-product between the velocity vector and the normal vector to the cross-sectional area was integrated throughout the volumetric ROI.
A workflow chart with the description steps of the 4D ultrafast ultrasound flow imaging is given in figure 2.

Experimental methods
2.2.1. 4D ultrafast imaging hardware. A customized, programmable, 1024-channel ultrasound system (Provost et al 2014) was used to drive a 32-by-32 matrix-array probe centered at 8 MHz with a 90% bandwidth at −6 dB, a 0.3 mm pitch and a 0.3 mm element size (Vermon, France). Briefly, the system was composed of four Aixplorer systems (Supersonic Imagine, France), each of which provided 256 transmit channels and 128 multiplexed receive channels. The four systems were assembled and synchronized which resulted in 1024 channels in transmission and 512 multiplexed channels in receive. Since the receive channels were multiplexed to 1 of 2 transducers elements, each emission was repeated twice to synthetize 1024 channels in reception. Specifically, data from two consecutive identical emissions were concatenated to synthetically create one RF data set, which corresponded approximately to an acquisition in which all transducers had received simultaneously, effectively dividing the maximum PRF by 2.
2.2.2. Imaging sequences parameters. The imaging sequences used for both the in vitro and the in vivo studies were built with similar transmit/receive parameters. Each plane wave was defined by a pair of angles between the normal vector of the plane wave and the x-and y-directions of the matrix-array probe. In this study, all the sequences used the following pairs in transmission: (−7°,0°) and (0°,7°) at 5 MHz frequency. The number of emissions was small to achieve a high volume-rate (i.e. the PRF divided by N) to limit aliasing effects. The angles used were chosen to be as large as possible while achieving a 8 × 8 mm 2 field-of-view (in which both waves propagate) at the 25 mm depth. The following nine receive angles: were used. The f/D ratio was chosen constant throughout the volume, and was set to its smallest possible value in this geometry, i.e. 5. The effective PRF was 8050 Hz, which corresponded to a volume rate equal to 4030 Hz. The acquisitions lasted 124 ms in vitro and 1.24 s in vivo, which corresponded to 500 and 5000 volumes, respectively.

In vitro experiments: a vessel phantom study.
Experiments to evaluate 3D velocities estimation in a full volume were conducted using a silicon tube-vessel phantom with a 4 mm internal diameter embedded at a 20 mm depth in an 1% agar, 99% propanol water-based gel. A gear pump (model 700-D, ATS Laboratories Inc., USA) supplied steady-state flow circulation of blood-mimicking fluid (model 707, ATS Laboratories Inc., USA) at flow rates controlled by an industrial flowmeter grade (model 700-D, ATS Laboratories Inc., USA), with a 2% reading accuracy. An 82° angle between the probe surface and the tube was measured on the B-mode volume. The technique was evaluated for six different volumetric flow rates. Table 1 indicates the volumetric flow rates used, along with the corresponding maximum peak velocities, but also the Nyquist limit criteria.

In vivo experiments: common carotid artery and carotid artery bifurcation
assessment. The in vivo feasibility was assessed by imaging the human carotid artery in two healthy volunteers. The ultrasonic sequences complied with the Food and Drug Administration (FDA) requirements (510k Track 3, FDA) regarding the mechanical index (MI) and the spatial-peak time average (I SPTA ) with a derating factor of 0.3 dB cm −1 MHz −1 . The MI 0.3 and I SPTA0.3 were 0.2 and 96.8 mW cm −2 , respectively, when assuming a repetition time of 1 s. The temperature increase at the surface of the probe was inferior to temperature uncertainty of the measurement system (±0.1 °C) during the test sequences.
The experiments were conducted and performed by a trained medical doctor. Prior to the 4D ultrafast acquisitions, a real time bi-plane B-mode imaging was used to achieve correct positioning and visualization of the carotid artery vasculature. Simultaneous ECG acquisitions were also performed, and all results were presented with a time origin corresponding to the R-wave.
For the common carotid artery measurements, both the 4D velocity vector distribution and the volumetric flow rate were evaluated in both subjects.
2.2.5. Data processing and analysis. All beamforming and velocity vector estimation processing were performed offline using a graphical processing unit (Tesla K40, NVidia, USA) within a Matlab (The MathWorks Inc., USA) interface. The in vitro and in vivo acquisitions required a total post-processing time of approximately 15 min and 1 h, respectively. This duration included all 4D ultrafast data transfers and calculations, i.e. the beamforming, clutter filtering, velocity vector components estimation, segmentation, and volumetric flow rate calcul ation. Calculation times were dependent on imaging depth, imaging sampling, number of transmit/receive angles and the total number of acquired volumes.
A trained medical doctor performed the segmentation used for the calculation of the volumetric flow rate on B-mode volumes reconstructed using a conventional coherent compounding delay-and-sum beamforming (Provost et al 2014) applied onto the same data used for 4D ultrafast flow imaging. 4D rendering of the 4D ultrafast B-mode volumes and 4D smoothed velocity vector maps were performed using the volume rendering (i.e. 'volren') and the illuminated streamlines functions of the Amira software (6.0.1, Visualization Sciences group, MA, USA), respectively.
Finally, in order to evaluate the precision of the velocity vector estimation, a normalized autocorrelation coefficient r was computed for each transmit-receive pairs, as follows: where Z is the axial voxel kernel-size in the z-direction (Z was equal to 5 in this study). Figure 3 presents the results of the in vitro experiment. Figure 3(a) shows the experimental apparatus. Figure 3(b) shows the velocity vector field rendered using streamlines overlaid onto the B-mode volumes for one representative volumetric flow rate (160 ml min −1 ). The streamlines' orientation indicated a flow direction aligned with the tube axis and their color showed velocities that are higher in the center and decrease toward the edges of the vessel, which is consistent with a parabolic flow pattern. Indeed, parabolic profiles can be observed in figure 3(c) for two representative volumetric flow rates (100 and 290 ml min −1 ). Note that parabolic patterns are observed in the surface velocity profiles, but also in the velocity profiles obtained in the center of the tube where larger velocities are observed along the y-direction, indicating the main direction of the flow. Table 1 shows the volumetric flow rate estimation, their respective mean estimated-error and mean auto-correlation coefficients for all the measurements.

In vitro experiment: vessel-phantom validation
The volumetric flow rates estimation was in good agreement with the ones measured experimentally by an industrial flowmeter, from 100 to 290 ml min −1 . At 360 ml min −1 , the averaged volumetric flow rate estimation error was approximately 5%. At 490 ml min −1 , the error made on the flow rate reached 18.5%, which is expected since the associated maximum peak velocity of 1.3 m s −1 exceeded the Nyquist velocity (i.e. 1.2 m s −1 ) by 8.4%. These measurements were associated with a decrease of the average auto-correlation coefficient to 0.70.   Figure 4(c) shows the average velocity inside each of the three ROIs of size (1.1 × 1.1 × 0.5) mm 3 . ROI A, B, and C were located near the anterior wall, in the center, and near the posterior wall of the artery, respectively.

In vivo experiment: common carotid artery and carotid artery bifurcation assessment
Movie 1 (available at: stacks.iop.org/PMB/61/L48/mmedia) shows the corresponding 3D velocity distribution during an entire cardiac cycle. The direction of the flow was aligned with the vessel axis and the variations observed in the velocity magnitude were consistent with the phases of the cardiac cycle.     Figure 5 presents the results of the in vivo carotid artery bifurcation experiments. Figure 5(a) shows the positioning of the ultrasound probe used to image the carotid artery bifurcation. Figure 5(b) shows the velocity vector maps at two time points, i.e. at t = 0.4 s and t = 1.1 s and One can note that while a laminar flow is observed at 0.4 s, a vortex flow appeared at 1.1 s, i.e. immediately after peak systole, in the posterior part of the internal carotid artery (ICA). Figure 5(c) shows the average velocity inside each of the three ROIs of size (1.1 × 1.1 × 0.5) mm 3 . ROI A, B, and C were located in the external carotid artery (ECA), in the center of the bifurcation and exist of the CCA, and in the ICA, respectively. It should be noted that the velocity profiles varied from one ROI to another due to the different sizes and orientations of the branches, which in turn lead to different volumetric flow rates.
Movie 2 (available at: stacks.iop.org/PMB/61/L48/mmedia) shows the corresponding 3D velocity distribution during an entire cardiac cycle. The flow is directed toward the external and ICA branches. Flow patterns differed during the cardiac cycle. Indeed, while laminar flow occurred during diastole, a vortex was visible during systole in the ICA branch.

Discussion
In this study, we introduced and demonstrated the feasibility of performing 4D ultrafast ultrasound flow imaging. More specifically, we have developed the experimental methods for volumetric blood flow imaging, evaluated the technique accuracy and precision in vitro in a tube-vessel phantom, and we have showed the in vivo feasibility in the human carotid artery.
We have imaged the complete volumetric blood velocity vector distribution using 2D steered plane wave emissions at high frame-rates and a 3D multi-angle cross-beam Doppler method, in which a large number of parameter combinations was possible. Indeed, the number and orientation of the transmit/receive angles, as well as the f/D ratio could be chosen. We used a sequence of two 7° plane wave emissions. The small number of emissions (i.e. two) allowed for the imaging without aliasing of the expected maximum peak velocities distribution of the targeted clinical application, i.e. the carotid artery. The tilt angle (i.e. 7°) was chosen to be as large as possible, within element directivity limits, to enhance the precision of the velocity vector (Yiu et al 2014) while allowing for a sufficiently large field of view (i.e. 8 × 8 mm 2 ) in which all waves propagated at the required 25 mm depth. We used nine angles with a f/D ratio equal to 5 as a compromise between the sensitivity, which is linked to the subaperture size, the precision, linked to the tilt angle, and the computational considerations, i.e. memory and speed. A complete and quantitative optimization of this large parameter space is beyond the scope of this study but is the subject of on-going work. 4D ultrafast ultrasound flow imaging was evaluated in vitro in a silicon tube-vessel phantom at different volumetric flow rates. An excellent agreement with the theoretical velocity profiles and volumetric flow rates was found in the axial cross-section for velocities below the Nyquist velocity. A flow rate error of less than approximately 5% was achieved for flow rates below 360 ml min −1 , which corresponded to a velocity of 100 cm s −1 . When exceeding the Nyquist velocity, e.g. when the flow peak velocity exceeded 1.2 m s −1 , the error increased rapidly to 18.3%, which was also associated to a reduced auto-correlation coefficient.
4D ultrafast ultrasound flow imaging was then performed in vivo in the carotid artery of two healthy volunteers. The peak velocities were found to be in agreement with the velocities measured in previous studies using magnetic resonance imaging (MRI) (Harloff et al 2009) or with 2D ultrasound-based techniques (Holdsworth et al 1999, Udesen et al 2007, Jensen 2014, Fadnes et al 2015. Similarly, the volumetric flow rates obtained were found to be in good agreement with the ones reported in the literature (Steinman et al 2002, Marshall et al 2004, Likittanasombut et al 2006. Finally, we imaged turbulent flow patterns at the carotid artery bifurcation. While it has been shown by 2D ultrasound velocity vector imaging techniques that such vortices can appear in the ICA during approximately 100 ms before the onset of diastole (Udesen et al 2007), it had not been imaged in 4 dimensions using ultrasound imaging.
4D ultrafast blood flow imaging may allow for the operator-independent, non-invasive, and direct determination of arterial volumetric flow rate. In clinical practice, the estimation of volumetric flow rates cannot be performed accurately with ultrasound imaging and is currently estimated by indirect measurements including the luminal cross-sectional area and the mean spatial velocity that induce important estimation errors (Hoskins 1990). Moreover, in contrast to MR flow imaging, 4D ultrafast ultrasound flow imaging could be used in a daily practice and at bedside. 4D ultrafast ultrasound flow imaging was developed for the angle-independent mapping of the 3D blood flow velocities distribution in contrast to the clinically used conventional angledependent ultrasound Doppler imaging. In a recent study, in vitro 3D velocity vector mapping was achieved with high temporal resolution by a 3D transverse-oscillations technique , but it was limited to two 2D planes rather than obtained in an entire 3D volume. Also, the V-echo-PIV technique (Falahatpisheh and Kheradvar 2016) provided 3D maps of all three components of the velocity vector in the heart in vivo, but this technique was limited by the imaging volume rate (i.e. 30 volumes s −1 ) and a contrast-agents injection requisite. In order to achieve 4D blood flow imaging, a 4D ultrafast ultrasound imaging prototype with 1024 electronic channels was used in this study. Such a high channel count system is required to perform 2D tilted plane wave emissions at high frame rates and 3D multi-angle crossbeam beamforming. To the best of our knowledge, only two research systems have cur rently these capabilities (Jensen et al 2013, Provost et al 2014. We can envision however, that such systems will be increasingly available with the recent development of high channel count electronics commercially available (Flynn et al 2011).
Current limitations include the choice and the number of the transmitted tilted angles to the specific clinical application and imaging field-of-view, and the need for manual segmentation for the determination of the volumetric flow rate. However, automatic segmentation algorithms could be used to provide a fully user-independent flow rate measurement. In order to improve the precision of high velocity estimation, the volumetric rate would need to be increased, which could be done using 1024 simultaneous electronic receive channels, which would double the Nyquist velocity. Inversely, it would also allow the use of a larger number of plane wave emissions to increase the precision of the velocity vector estimation.

Conclusion
4D ultrafast ultrasound flow imaging can map in 3D velocity vectors of blood flow in large arteries at high temporal resolution. Its precision and accuracy were respectively evaluated both in terms of velocity and volumetric flow rate measurements in vitro and its feasibility was demonstrated in vivo in the human carotid artery. 4D ultrafast ultrasound flow imaging could become a diagnostic imaging tool for the real-time evaluation of arterial blood flows in clinical practice.