SPATIO-TEMPORAL FLOW AND WALL SHEAR STRESS MAPPING BASED ON INCOHERENT ENSEMBLE-CORRELATION OF ULTRAFAST CONTRAST ENHANCED ULTRASOUND IMAGES

In this study, a technique for high-frame-rate ultrasound imaging velocimetry (UIV) is extended first to provide more robust quantitative flow velocity mapping using ensemble correlation of images without coherent compounding, and second to generate spatio-temporal wall shear stress (WSS) distribution. A simulation model, which couples the ultrasound simulator with analytical flow solution, was implemented to evaluate its accuracy. It is shown that the proposed approach can reduce errors in velocity estimation by up to 10-fold in comparison with the coherent correlation approach. Mean errors (ME) of 3.2% and 8.6% were estimated under a steady flow condition, while 3.0% and 10.6% were found under a pulsatile condition for the velocity and wall shear rate (WSR) measurement, respectively. Appropriate filter parameters were selected to constrain the velocity profiles before WSR estimations and the effects of incorrect wall tracking were quantified under a controlled environment. Although accurate wall tracking is found to be critical in WSR measurement (as a 200 μm deviation from the wall may yield up to a 60% error), this can be mitigated by HFR imaging (of up to 10 kHz) with contrast agents, which allow for improved differentiation of the wall-fluid boundaries. In vitro investigations on two carotid bifurcation phantoms, normal and diseased, were conducted, and their relative differences in terms of the flow patterns and WSR distribution were demonstrated. It is shown that high-frame-rate UIV technique can be a non-invasive tool to measure quantitatively the spatio-temporal velocity and WSS distribution. (E-mail: mengxing.tang@imperial.ac.uk) © 2017 The Author(s). Published by Elsevier Inc. on behalf of World Federation for Ultrasound in Medicine & Biology. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


INTRODUCTION
Studies have shown that hemodynamic shear stresses have a strong influence on the initiation and development of various vascular diseases. In response to imposed mechanical forces, endothelial cell layers not only endure the morphologic and functional changes (Chistiakov et al. 2017;Chiu and Chien 2011), but also trigger the biochemical event that results in the deposition and adherence of the platelet or the formation of atheroma (Caro 2009;Cecchi et al. 2011;Dhawan et al. 2010). Such observations are convincing evidence that local hemodynamic factors are associated with vascular diseases. It is therefore of clinical significance to quantitatively measure the flow patterns and wall shear stress (WSS), using the existing medical modalities.
Blood flow patterns, flow velocities or flow volume have been widely used to describe the spatio-temporal hemodynamic within a particular vessel or lesion. Although several techniques of measuring the flow velocities exist clinically, a standard technical basis for fluid shear estimation in vivo has not been established (Katritsis et al. 2007;Papaioannou and Stefanadis 2005). The most convenient and straightforward way to estimate WSS is based on the Hagen-Poiseuille formula, where the shear stress in a relatively straight vessel is computed in terms of either the average or maximum velocity and the radius of the vessel lumen (Mynard et al. 2013). However, the accuracy of such expression is dependent on the validity of the following assumptions: (i) The flow is steady and of a parabolic in nature; (ii) the blood behaves as a Newtonian fluid with constant viscosity; (iii) the vessel is axisymmetric; and (iv) the vessel wall is rigid. These conditions are met minimally in the physiologic environment where the flow is mainly pulsatile, the vessel is distensible and its diameter is inconsistent. These conditions therefore limit the application of such a technique to estimate the WSS. Another technique to estimate the WSS is based on the determination of the wall shear rate (WSR), which is the flow velocity gradient near the vessel wall (Cheng et al. 2002;Poelma et al. 2012). This produces a more reliable WSR measurement with no assumption of the geometry and flow required; however, the accuracy of such a technique is dependent on the temporal and spatial resolution of the modality used to obtain the flow velocity profile and track the wall position, and the interpolation or extrapolation algorithm. Because of such constraints, in vivo WSS measurement is limited and relies primarily on computational fluid dynamics (CFD). In this context, numerical simulation is employed to investigate the WSS distribution of a given geometry measured by means of imaging modalities. However, the accuracy of the simulation can be affected by the underlying assumptions on the geometry, wall properties, fluid properties and most importantly the initial and boundary conditions.
To calculate the WSS in vivo, flow velocity information can be assessed using existing non-invasive imaging modalities such as phase contrast magnetic resonance (Masaryk et al. 1999;Stalder et al. 2008) and ultrasound (Poelma et al. 2012;Zhang et al. 2011). However, it has been a challenge clinically to provide a reliable WSS distribution with the existing medical imaging technique.
A high-frame-rate (HFR) plane wave ultrasound imaging velocimetry (UIV) system has been developed and evaluated recently for the generation of flow velocity mapping with high spatial and temporal resolution (Leow et al. 2015a). This technique benefits from the combination of the HFR ultrasound imaging, contrast imaging and UIV algorithm. HFR ultrasound imaging provides high temporal resolution, of up to tens of thousands of frames per second, to track high spatio-temporal variation flow within the field of view. The use of microbubbles-together with a pulse-inversion contrast-imaging sequence-enhances the signal-to-noise ratio (SNR) and helps differentiate between the wall-fluid boundaries. The high temporal resolution and SNR imaging technique therefore allows the UIV algorithm to estimate the flow velocities and subsequently the WSS unambiguously.
Here, a previous study has been extended to provide both robust quantitative velocity flow field measurement and spatio-temporal WSS, using a modified UIV technique based on ensemble-correlation of images without coherent compounding. An incoherent ensemble-correlation approach such as this avoids the effects of motion on the coherently compounded images. The accuracy of the system in estimating flow velocity vectors and WSR mapping is evaluated. Ultrasound flow simulation, which consists of a straight tube phantom driven by steady and pulsatile flow with known reference values, fully evaluates this technique. In vitro studies on two anthropomorphic carotid bifurcation phantoms are also conducted to further highlight the performance of the technique.

METHODS AND MATERIALS
Ultrasound flow simulation Flow model generation. Synthetic ultrasound images were generated from ultrasound flow simulation that coupled the analytical flow solution with the ultrasound simulator. Such realistic models are necessary to compare the developed technique with the known reference value computed from the analytical solution. Under such circumstances, a rigid model, consisting of a straight cylindrical vessel phantom filled with Newtonian fluids driven by an analytical solution, was created.

Poiseuille flow
Straight vessel phantoms with inner diameters of 3 mm and 6 mm were created. Established steady laminar flow with a parabolic flow profile was generated by moving the flow scatterers at radius r from the center of the vessel with the velocity as follows: where v0 is the centreline velocity and R is the radius of the vessel. The WSS, given a fluid with dynamic viscosity µ, is defined as: In this study, flow simulations under three flow rates (v0 = 20, 50, 80 cm/s) were generated in a 3-mm diameter vessel to investigate the accuracy of flow estimation. The vessel was positioned with a beam-to-flow angle of 60°and the radial velocity errors were calculated. In a 6-mm diameter vessel, steady flow with the centerline velocity of 50 cm/s, Reynolds number (Re) of 1500 and resulting WSR of 333 s −1 was also simulated. The vessel was oriented with a beam-to-flow angle of 90°for the simplicity of the WSS analysis.

Womersley flow
While arterial flow is pulsatile, Womersley flow was generated as described previously (Evans and McDicken 2000). In brief, pulsatile flow is considered as the sum of both steady flow component and a series of oscillatory components. Poiseulle's equation can be reformulated to describe the relationship between pressure difference and volume flow rate (Womersley 1955): where ΔP t sin ω φ + ( ) is the pressure difference, α is the Womersley number and M′ 10 and ε′ 10 are complex functions of the non-dimensional parameter α ω υ = R . From the volume flow rate, it is possible to calculate the velocity profile from a sinusoidal flow (Evans 1982 where Jm is the m th order Bessel function of the first kind, τ α = j 3 2 , ψn and χn represent the amplitude and angle of the complex function ψn. Therefore, the velocity profile can be calculated from volume flow rate when the entrance effects are neglected.
Following the Fourier decomposition theorem, the mean velocity waveform can be expressed in a series of sinusoidal waveform such that V t V V cos n t n n n The time evolution of the velocity profile can be determined by summing the velocity profile from each harmonic.
In this study, the velocity profile from the common carotid artery (CCA) and the common femoral artery (CFA) were generated using the data presented in Table 1.

Ultrasound simulation
Field II (Jensen 1996;Jensen and Svendsen 1992) was used to simulate a realistic ultrasound system. Arbitrary transducers and imaging pulse sequence can be created on the basis of the linear systems theory and the spatial impulse responses estimation described previously (Stepanishen 1971;Tupholme 1969). In this context, an ultrasound field at any given point can be estimated by taking into account the following: transducer estimation function, spatial impulse response of the transmitting and receiving aperture and the electro-mechanical transfer function of the transducer.
In Field II, tissues and blood are modelled as a collection of point scatterers with various reflectivity. Ten scatterers per resolution cell were simulated to ensure a fully developed speckle with Gaussian distributed radio frequency (RF) signal (Thijssen 2003). The rigid vessel walls were generated as stationary scatterers while the Newtonian fluids were represented as random distributed point scatterers whose positions were updated for each transmitted ultrasound beam during the simulation, using the displacement determined by the analytical solution.
The parameters of a linear array transducer (L12-3 v, Verasonic, Redmond, WA, USA) and compounded plane wave imaging sequence were used to simulate the RF  signal, as presented in Table 2, to represent the ultrafast ultrasound system used in the experimental study. A higher sampling frequency is necessary to reduce the side lobes and better represent the simulated spatial impulse response. The transducer elements were also divided into smaller sub-apertures such that every scatterer is in the far field of each sub-element (Jensen and Svendsen 1992). In this case, the elements were further divided into four subelements in the elevational direction to ensure the far field of each sub-element was at an axial distance of 1.3 mm while keeping the simulation time short. Also to simulate a realistic ultrasound field, random Gaussian noise was included in individual pre-beamformed RF-data-frame yielding an SNR of 20 dB. A software graphics processing unit (GPU) beamformer (Yiu et al. 2011), in which delay and sum operations were applied to the RF data, was used to generate synthetic ultrasound images frame for further analysis.

In vitro flow experiment
Ultrafast plane wave imaging system. Aresearch platform (Vantage 128, Verasonic,) was configured to acquire HFR plane wave ultrasound images.Alinear array transducer (L12-3 v, Verasonic) was mounted to image the flow phantom. Real time RF data were recorded in the local memory and later transferred to a computer through a highspeed Peripheral Component Interconnect Express (PCI-Express) connection and beamformed using the GPU beamformer for further analysis. This system has been used in previous flow velocity mapping studies (Leow et al. 2015a).
Microbubble contrast agents and contrast imaging mode. Microbubble contrast agents have been used in this study to enhance the ultrasound signal within the flow stream. It was prepared using a recipe previously described in Sheeran et al. (2011), resulting in a solution contained of 5 × 10 9 micobubbles/mL with an average size of 1 µm (Sennoga et al. 2012). In this study, microbubbles were diluted in blood-mimicking fluid (BMF), which was composed of 90% pure water and 10% glycerol, to a concentration of 2 × 10 5 microbubbles/mL. This concentration is clinically relevant and has been used in previous studies (Leow et al. 2015a(Leow et al. , 2015bTang et al. 2010).
Pulse inversion (PI) imaging sequences were implemented to acquire ultrasound contrast images. Although it reduces the imaging frame rate by the number of pulses transmitted to form an image, the developed system is capable of keeping the temporal resolution in the kHz range when coupled with HFR imaging. In addition, the PI images are found to be superior to either B-mode or second harmonic imaging in terms of its specificity in imaging microbubbles (Stanziola et al. 2016). It is capable of suppressing tissue signals that oscillate primarily at fundamental frequency, while maintaining harmonic signals generated from microbubbles. PI imaging therefore enhances the contrast and SNR of the generated image.
In this study, compounded plane wave PI imaging sequences were implemented using the Verasonics research platform. Six plane waves with angles spanning −10°to 10°were transmitted at the 10-kHz pulse repetition frequency. For each angle, the L12-3 v transducer (Verasonics) was driven to transmit a 4 MHz 1-cycle pulse followed by its phase-inverted counterpart before transmittin from the next angle. The RF echoes were then transferred back to a computer for further analysis. Typical settings for the data acquisition are presented in Table 3.

Velocity estimation using a modified-UIV analysis.
Advanced UIV analysis with an iterative window deformation model (Scarano 2002) was performed to estimate the flow velocity. Such a technique is implemented in a hierarchical scheme that progressively compensates the inplane flow motion on a smaller scale. In brief, crosscorrelations were performed iteratively with progressive grid refinement to estimate the local flow field. For each iteration, a three-point Gaussian peak fitting and universal outlier detection (Westerweel and Scarano 2005) were implemented to estimate sub-pixel displacement and eliminate spurious vectors. The resulting estimations from the previous iteration are interpolated to build a displacement predictor that will be used to deform the images for the next iteration until a convergence criterion is satisfied. The validated velocity field with a higher resolution is displayed as vector velocity mapping. An ensemble correlation approach, in which multiple correlation functions computed from several image pairs are averaged before peak searching in the correlation plane, can be implemented. This serves to further enhance the tracking accuracy. Such an approach has been widely used in the optical field to achieve velocity tracking at the single pixel level (Wereley and Meinhart 2010), also to increase the SNR of a correlation function to reduce correlation errors (Meinhart et al. 2000;Poelma et al. 2011).
In this study, an incoherent ensemble correlation approach is developed. A joint correlation function is obtained by summation of the cross-correlation maps acquired from the multiple low-resolution images. This is different from the coherent ensemble correlation approach, which is done on the final coherent compounded images. Such a modification is performed as the decorrelation of the speckle pattern depends not only on the velocity gradient, but also on the motion artefacts that generated from the flow motion between pulse emissions (Denarie et al. 2013). By taking the average of the correlation plane computed from the low-resolutions images, the motion artefacts can be circumvented and robust velocity estimation can be achieved.
The schematic workflow of modified-UIV analysis and the principle of incoherent ensemble-correlation are illustrated in Figure 1. Comparisons among the incoherent ensemble-correlation approaches with the coherent ensemble correlation method were demonstrated in the in vitro experiment to evaluate the new method.
WSS estimation. WSS can be measured directly from the velocity profile. This allows the quantification of WSR not only in a simple straight tube geometry, but also in complex configurations such as a bent vessel and a bifurcation. In this study, WSR is assessed by estimating the velocity gradient from the velocity profile tangential to the wall. It should be noted that because shear rate is a gradient, it can be easily affected by UIV estimation errors and noises. A Savitzky-Golay (SG) filter is introduced locally to smooth the velocity profile to mitigate such effects before shear rate estimation. This is a digital filter that smooths the data by fitting a subset of adjacent data points with a low-degree polynomial by linear least-squares methods (Savitzky and Golay 1964). Prudent steps need to be taken to avoid underfitting or overfitting the flow profile. Hence, the effects of the filter parameters chosen, in terms of the polynomial order and the length of the convolution function, on the accuracy of WSR estimation were investigated.
While a non-slip boundary condition was applied to estimate the velocity gradient, accurate identification of the vessel wall is essential. In the experimental study, a localized region-based active contour segmentation algorithm (Lankton and Tannenbaum 2008) was implemented to dynamically track the vessel wall before the velocity estimation. It is a robust and accurate technique that uses local image statistics to evolve a contour for object segmentation. However, such segmentation is not implemented to the steady flow simulation study. The wall positions are deliberately deviated from the reference wall position in the steady flow simulation study to investigate the effect of segmentation errors toward the WSR accuracy.
By assuming the working fluid is a Newtonian fluid, WSS can be determined from the shear rate tensor as follows: where µ is the dynamic viscosity and ε′ 12 is the tangential component shear rate tensor. To account for any misalignment between the image coordinate and the coordinate along the wall, the shear rate tensor was computed by transforming the original strain rate tensor in the 2-D Cartesian coordinate to the wall-oriented coordinate system as shown in Figure 2 (Charonko et al. 2009;Kundu 2012).
The transformation can be written as where εij is the original strain tensor in the image Cartesian coordinate defined as Cij is the 2-D rotation transformation matrix, defined for each point along the wall in terms of rotation angle θ between the original and wall-oriented coordinate, as follows C cos sin

Error quantification
The purpose of this study is to evaluate the performance of the developed technique, and two error metricsthe mean error (ME) and the normalized root mean square error (NRMSE)-were considered, as described in eqns (12) and (13). where y and yt are the UIV estimated value, the reference value, y t max ( ) is the maximum reference value and n is the number of time sample.
Anthropomorphic flow phantoms. Two anatomically realistic flow phantoms (Lai et al. 2013), normal bifurcation geometry and diseased bifurcation geometry with 50% eccentric stenosis were investigated to assess the practicality and feasibility of this technique to provide both velocity mapping and WSR. The phantoms contain wall-less vessels and are fabricated with polyvinyl alcohol cryogel (PVA-C) as the tissue-mimicking material. They contain three vessel branches, which mimic a carotid bifurcation junction in the vicinity of CCA, internal carotid artery (ICA) and external carotid artery (ECA). It is noted that the inner diameters of CCA, ICA and ECA are 6, 4.2 and 3.2 mm, respectively.
The inlet of the carotid model was connected to a pulsatile pump (1405, Harvard Apparatus, Holliston, MA, USA). The circulating fluid comprises BMF solution and diluted microbubbles. It is circulated at a pulse rate of 60 pulses/min and an average flow rate of 2 mL/s with a 4 mL/s systolic flow rate to reproduce physiologic conditions.

Ultrasound flow simulation
Steady flow. The coherent and incoherent ensemblecorrelation approaches of the UIV algorithm were performed to estimate the 2-D flow velocities as illustrated in Figure 3. The radial velocity estimation errors are summarized in Table 4. It is observed that under slow flow conditions (v0 = 20 cm/s), an incoherent ensemble correlation approach performs as well as the coherent ensemble correlation approach. However, when the flow velocity increases, the accuracy of the coherent approach degrades and significant errors are quantified when the flow is fast (v0 = 80 cm/s). The incoherent ensemble correlation approach, on the other hand, is more robust against the motion artefacts and highly accurate measurement is estimated under all simulated flow conditions.
With the known vessel wall location, simulated images were segmented and subsequently estimated using the UIV algorithm. The 2-D velocity mapping, averaged velocity profile and shear rate profile are illustrated in Figure 4. Consistent parabolic profiles, as represented by color-coded velocity vectors, can be seen across the straight vessel. An ME of −3.2% or NRMSE of 1.9% was found when comparing the estimated velocity profile to the analytical velocity profile. Despite the high accuracy of the velocity measurement, a small fluctuation in the velocity profile may affect the WSR estimation as shear rate is estimated from the velocity gradient, as shown in Figures 4b and 4c. The shear stress error seems to be greater at the -r position (Fig. 4c) because of noise near the boundary. The velocity is close to zero near the boundary; therefore, the same amount of noise results in more error in the gradient measurement computed directly from the pre-filter velocity profile. However, with an SG filter applied to smooth the velocity profile before the shear rate estimation, a more robust WSR measurement is obtained.
The effects of the filter length and the polynomial order were investigated to optimize the filter settings, as shown in Figure 4d. It should be noted that only second and third order filters were investigated because a highorder polynomial is generally less robust and sensitive to small-amplitude and random distributed error (Lou et al. 1993). The findings, as shown in Figure 4d, suggest that an increase in the filter length relative to the diameter of the vessel produced less error with a smaller variance; however, it may intrinsically restrict the flow profile to certain polynomial shape depending on the polynomial order selected. This forces the velocity profile to comply more with the velocity profile in the region farther away from the wall than to the profile near the wall. Hence, filter parameters were chosen such that the mean estimation error and variance are less than 10%. In the later study, we implemented a second order 0.2-D filter, resulting in a 6.9 ± 5.8% error (D = vessel radial diameter), and a third order 0.3-D filter, resulting in an 8.6 ± 6.8% error.
Using the optimized filter parameters, the WSR estimation error in a function of estimated wall position relative to the reference wall position is illustrated in Figure 5a. Regardless of the filter order, wall shear rates are overestimated when the wall positions were incorrectly localized within the vessel and underestimated when the wall position were localized outer from the reference wall position. Such observation can be explained by Figure 5b, which plotted the shear rate profile corresponding to various wall locations identified. The wall shear rates, estimated at both ends of the shear rate profile, diverged from the reference profile up to 60% when an inaccurate wall position of 200 µm was identified. Such discrepancies can be accounted for the by the low-pass filtering effect of the correlation windows (Willert and Gharib 1991) and the non-slip boundary condition, which imposes a zero velocity at the boundary.
Pulsatile flow. The high spatio-temporal flow velocity vectors and WSR mapping of the CFA and CCA are illustrated in Supplementary Video S1 and S2 (online only, available at https://doi.org/10.1016/j.ultrasmedbio .2017.08.930). Both videos are generated with a simulated acquisition frame rate of 1 kHz but played back at 100 Hz. To highlight the flow structures in CCA and CFA, Figures 6 and 7 show the temporal snapshot of the velocity vectors and WSR mapping at various phases of the pulsatile flow indicated in Figures 6c and 7i. Spatial average Fig. 2. The tensor transformation from the image-oriented Cartesian coordinate (x2,x2) to the wall-oriented coordinate system. Note that u is a 2-D vector, εij are the second-order tensors and ε′ 12 is the wall shear rate on the image-oriented coordinate system. velocity profiles overlaid on the analytical solutions are also demonstrated next to the temporal snapshots. Unidirectional flows, mainly parabolic, can be observed throughout a pulse cycle in CCA (Figs. 6a and 6b), while more complex flow structures, including non-parabolic flow and reverse flow, can be observed in CFA. The transition from an undeveloped flow with blunt profile (Fig. 7a) to a parabolic flow (Fig. 7b) can be seen during the systole  phase while a flow reversal (Fig. 7c), and a non-uniform flow profile (Fig. 7d) can be observed during the diastolic phase of the cardiac cycle.
Such high dynamic flow also created a temporal variation in the WSR overlaid on the vessel wall. The magnitude of the WSR varied according to near-wall flow  showing a strong agreement between the measurements and the references. From the quantitative measurement, the estimated centerline velocities were found to match the analytical solutions with ME or NRMSE of 3.0% or 2.8% in CCA and 1.1% or 0.3% in CFA. The spatialaveraged WSR are illustrated in Figure 8 to quantify the WSR estimation error. Only the upper WSRs are demonstrated because the WSR is symmetric. When comparing to the analytical solution, significant error can be observed when no SG filter is applied before estimating the WSR in the CCA and the CFA. Such an error is similar  to those observed in Figure 4 and can be accounted for by the variance in the velocity profile close to the boundary. The WSR estimation is significantly improved, as indicated in Figure 8, after filtering the velocity profile with a third order SG filter. The estimated spatial averaged WSR shows reliable agreement with the analytical solution and the estimation error were summarized in Table 5. These results also indicate a robust WSR estimation as the addition noise, or an SNR of 20 dB, does not change the measurement compared with the ideal situation.

In vitro results
Using the optimized settings obtained from the simulation study, Supplemenatry Video S3 and S4 (online only, available at https://doi.org/10.1016/j.ultrasmedbio .2017.08.930) illustrate the time-resolved quantitative visualization of a high dynamic flow that emerged within a healthy and a diseased carotid bifurcation phantom under the pulsatile flow condition. Both videos are acquired at 1.67 kHz but played back at 100 fps to render the fast temporal evolution of the flow patterns and the WSR mapping.
To highlight the spatio-temporally varying flow patterns induced by the two geometry configurations, a few temporal snapshots of the quantitative UIV measurement acquired at four cardiac phases are shown in Figure 9. For the healthy carotid phantom, a unidirectional forward flow can be observed during the systolic upstroke and peak systole (Fig. 9a). High-velocity flow up to 0.7 ms -1 and WSR up to 2000 s -1 can be observed throughout the carotid bifurcation phantom except at the ICA bulb where a relatively low velocity and WSR can be seen. A vortex subsequently emerges near the ICA sinus immediately after the peak systole and slowly dissipates during the post-systolic phase (Fig. 9b), causing a flow reversal and a negative WSR near the carotid bulb. In addition, it is worth noting that at the inner wall of proximal ICA and ECA, where the branching is located, a relatively high WSR can be observed because of the forward streamline flow. When the dicrotic pulse produced a secondary upstroke, the flow moves forward again, as shown in Figure 9c. This is followed by the reappearance of the vortex near the carotid bulb during the post-dicrotic phase, as shown in Figure 9d. However, during the diastole phase, the flow moves relatively slowly and therefore lower WSR is expected.
On the other hand, a distinct flow pattern can be observed in the diseased carotid phantom, particularly at the ICA branch because of the eccentric stenosis at that region. As a similar inlet flow is applied, the flow in the ECA branch of the diseased phantom is relatively consistent with normal phantom; however, the formation of a highvelocity flow jet at the ICA stenosis throat can be been clearly during the systolic upstroke. This follows a path from the inner wall at the proximal ICA to the straight outer-wall region at the distal ICA, thus creating a region with higher WSR. The reorientation of the jet flow tangential to the wall and the reflection of the jet stream on incidence with the outer wall can also be observed. In addition to the jet formation, the formation of the two vortices, one at the carotid bulb and the other along the inner wall at the proximal ICA, are observed in Supplementary   Fig. 8. Spatial-averaged WSR estimation overlay on the analytical solution (left) in a CCA and (right) a CFA. WSR = wall shear rate; CCA = common carotid artery; CFA = common femoral artery. CCA = common carotid artery; CFA = common femoral artery; ME = mean error; NRMSE = normalized root mean square error. Fig. 9. Key quantitative measurements derived from UIV to highlight the spatio-temporal varying flow profile and WSR in (a-d) a healthy carotid bifurcation phantom; (f-i) a diseased carotid bifurcation phantom. (e-j) Centerline velocity extracted from the common carotid artery. The relative position of each flow velocity and WSR mapping is presented in the velocity plot. UIV = ultrasound imaging velocimetry; WSR = wall shear rate.
Video S4 (online only, available at https://doi.org/10.1016/ j.ultrasmedbio.2017.08.930) and highlighted in Fig. 9f-9i. While the jet stream from the stenosed ICA bulb toward the outer wall at the distal end and the vortex at the carotid bulb can be seen clearly throughout the cardiac cycle, the vortex at the inner wall of the distal ICA emerge and dissipates twice within a pulse cycle. The recirculation at the distal ICA only occurs after the jet stream becomes narrower and crosses over from the inner wall to the distal ICA. This forms a flow divider between two recirculation zones during the post-systole and post dicrotic, as seen in Figures 9g and 9i. This is not the case during the period before peak systole and peak dicrotic wave, as shown in Figures 9f and 9h, as the jet is more dispersed after traveling past the inner wall. It is noted that the wall regions where the vortices are located constantly exhibit a low WSR in contrast to the wall area where the jet travelled. For example, the stenosis throat, the inner wall of the proxi-mal ICA and the outer wall of the distal ICA are subjected to high WSR. The relative difference between the normal and diseased, in terms of flow velocity and WSR patterns, can be clearly seen in Figure 9. Such findings agree with those previously observed (Kefayati et al. 2014;Yiu et al. 2014). Figure 10 shows the vector flow estimation between the classic ensemble correlation and the incoherent ensemble correlation. This illustrates the robustness of the developed technique. It can be clearly seen that during the peak systolic, the highly accelerated jet flow at the ICA throat, which cannot be tracked with the coherent ensemble correlation (N coherent compounded images use to performed correlation analysis), can now be estimated using the proposed method (N × M low resolution images). It was evident that the modified UIV is robust against motion artefacts, and a more accurate estimation can be achieved when the incoherent ensemble correlation is utilized.

DISCUSSION
The work presented in this study is an extension of previous work on an HFR-UIV system (Leow et al. 2015a). By correlating HFR images without coherent compounding, it eliminates the effects of motion on such compounding and allows more accurate quantification of the dynamic arterial flow velocity field. Furthermore, methods for quantifying spatio-temporal distribution of WSS are also developed. An ultrasound flow simulation, which integrates the analytical flow solution and the Field II ultrasound simulator, has been implemented to assess this novel technique. Such evaluation allows not only the assessment of the system's accuracy, but also the optimization of the system parameters.
From the simulation studies, the capability of the plane wave UIV system to provide a sensitive, accurate and fullfield of view velocity measurement is demonstrated in Figures 4, 6 and 7. The spatio-temporal flow dynamics of a pulsatile flow in terms of centerline velocity and flow velocity profile within the field of view are well estimated within the error of less than 10%. This is only possible with the use of high temporal resolution of plane wave UIV. It allows the measurement of a wider range of velocities, as well as the reduction of errors in UIV estimation that may arise in conventional line-by-line sweeping (Zhou et al. 2013).
Despite the high-accuracy velocity mapping, accurate measurement of WSR is not an easy task. WSR is derived from the velocity gradient near the wall; therefore, slight changes in the velocity gradient or the wall position may affect the estimation. An SG filter has been applied to smooth the velocity profile before estimating the WSR to mitigate the former problem. Such a filter has shown to reduce effectively the uncertainties in WSR measurement, as shown in Figure 4. Note that a careful selection of the order of the filter and coefficient length is necessary to avoid under-or overfitting estimation, similar to results where interpolation or extrapolation methods are used (Fatemi and Rittgers 1994). Generally, the order of the filter defines the shape of the fitted curve; whereas, the length of the filter determines the smoothness of the curve. The filter length can be chosen such that mean and variance of the WSR error start to plateau to avoid oversmoothing the flow profile to certain polynomial shapes. The filter order, on the other hand, can be chosen if the shape of the flow profile is known. The 2nd order filter slightly outperformed the third order filter as indicated in Figure 4d because the parabolic flow profile is better represented with a second order polynomial. However, physiologic flow profile is unnecessary in a parabolic shape as shown in the CFA in Supplementary Video S2 (online only, available at https://doi.org/10.1016/j.ultrasmedbio .2017.08.930) and Figure 7. The third order SG filter is therefore more appropriate as it provides a greater flexibility in smoothing the flow profile without compromising the WSR estimation accuracy even when the flow profile is parabolic. This is evidenced by the fact that the third order polynomial SG filter yields an accurate and robust WSR measurement across the full field of view both spatially and temporally as demonstrated in Figure 6 and 7.
While consistent measurement can be realized with filtering, a more stringent requirement to identify the wall position is also necessary for accurate WSR estimation. The effect of inaccurate wall tracking has been illustrated in Figure 5, where significant under-or overestimation may occur, depending on the wall location used in the WSR calculation. Such a condition can be attributed to the limited spatial resolution in imaging modalities as discussed in the literature, where a direct WSR measurement technique is implemented (Kim et al. 2004;Masaryk et al. 1999). However, this novel technique, which combined the plane wave ultrasound imaging with contrast imaging, may well overcome such difficulties. Not only do the microbubble contrast agents enhance the signal within the blood pool and increase the ultrasound imaging sensitivity, they also enable a higher specificity between vessel wall and flow region when coupled with a contrastspecific imaging sequence as shown in Figure 11. From 10 sequential images, the microbubbles and tissue regions were manually selected and the contrast-to-tissue ratio,

CTR
S S microbubble Tissue = ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ 20 10 log , was calculated. S microbubble is the mean microbubbles signal and S Tissue is the mean tissue signal computed from the red and yellow boxes indicated in Figure 11. The averaged CTR computed from 10 images were 11.5 ± 0.2 dB for the B-mode images and 19.2 ± 0.3 dB for the PI images, achieving a 7.7-dB improvement in this case. Such an elevated CTR in the PI imaging therefore allows a better wall segmentation in comparison with the B-mode imaging while plane wave imaging retains a high temporal resolution for wall tracking and velocity measurement. The capabilities of plane wave UIV to quantify and visualize complex flow dynamics within a physiologicrelevant geometry are also illustrated in laboratory experiments. The transient behavior that results from the flow pulsatility and the flow patterns because of the influence of vessel geometry can be clearly visualized in Supplemenatry Video S3 and S4 (online only, available at https://doi.org/10. 1016/j.ultrasmedbio.2017.08.930). The flow, which is fast, varies spatially and temporally and cannot be tracked using conventional UIV or standard Doppler, can now be visualized clearly and quantified using high-frame rate UIV.
Spatio-temporal WSR mapping, which is overlaid on the video, is another technical merit of this innovative technique. The WSR distribution seen in this study matches relatively well qualitatively with those optical PIV measurements demonstrated on a transparent and rigid carotid bifurcation model (Kefayati et al. 2014). However, it is accepted that a lower WSS/WSR magnitude is quantified here. This discrepancy can be attributed to several factors such as the variation in the flow rates, variation in fluid viscosity, uncertainties in the wall location, out-ofplane motion and the compliance of the model. Despite this, the relative variations among the WSRs of a healthy and a diseased carotid geometry can be clearly observed in Figure 9.
HFR imaging is able to increase the maximum velocity that can be measured; however, as the flow under sonification moves during the transmit event, coherent plane wave compounding is no longer achieved and may cause significant image artefact. This is dependent on the flow velocity and will limit the estimation accuracy among imaging frames as shown in Figures 3, 10a and 10c. A more robust velocity estimation technique is demonstrated using the incoherent ensemble correlation approach to circumvent the speckle decorrelation because of the image compounding. High-speed jet flow, which cannot be tracked using the classical correlation method, can now be estimated using this innovative approach, as illustrated in Figures 10b and 10d.
Speckle decorrelation may also arise from the outof-plane motion. The flowing microbubbles can move out of the imaging field of view after several imaging pulses. This will result in significant speckle decorrelation that affects the tracking of microbubbles, using the ensemble correlation method. The plane wave imaging technique, which is capable of imaging at ultrafast frame rate, is less affected by the out-of-plane motion. A more robust mea-surement is achieved when estimation is done using incoherent ensemble correlation as it is not suffering from the intensity or phase variations due to coherent summing of the echoes in the presence of in-plane or out-of-plane motion. While this study focused on the contrast enhanced HFR ultrasound imaging, there is a similar effect of motion within the compounded image frames on noncontrast HFR ultrasound imaging. The coherent ensemblecorrelation approach, with an integrated motion correction technique to remove motion artefacts in the coherent compounded images, has been demonstrated in the speckletracking technique without contrast agents (Joos et al. 2016). Contrary to this, the incoherent ensemble-correlation approach we proposed is a new approach that can enhance the robustness of the motion estimation without motion correction. It circumvents the effects of image artefact and increasing the SNR of the correlation function by averaging the correlation function computed from multiple lowresolution images. While in principle it can be applied to B-mode (without contrast agents) images as well, the sensitivity without contrast agents will be lower and the accuracy of the estimation will likely degrade. More studies are needed to optimize the proposed technique for noncontrast flow quantification.
In recent years, many new Doppler techniques have been developed, along with the advancement of ultrafast ultrasound imaging to overcome the conventional Doppler limitations. Such techniques include multi-beamed vector Doppler and transverse oscillation (Jensen et al. 2016a(Jensen et al. , 2016b. Although promising results to track twodimensional flow and the ability for real-time estimation have been demonstrated, Doppler techniques are still limited by aliasing, which is inherent in the velocity estimator (Swillens et al. 2010). Our approach, on the other hand, Fig. 11. The comparison between B-mode imaging and PI imaging demonstrates a higher specificity in imaging microbubbles and CTR. Note that tissue is highlighted in yellow and microbubbles (depicted in red) boxes are selected to compute CTR. PI = pulse inversion; CTR = contrast-to-noise ratio.
can be a more consistent and practical approach that is not only capable of estimating a wider range of velocities, but also robust against the motion artefacts. This, however, comes with the cost of the relatively high computational cost compared with the Doppler flow estimation approach. The processing time, from the UIV tracking to the WSR estimation, was 5.8 s for each measurement. The majority of the computational cost can be attributed to the large-scale cross-correlation analysis and iterative image deformation in the UIV algorithm, which works on multiple low-resolution images. With the evolution in graphical processing unit processing, the computation speed of the proposed method can be significantly increased and realtime processing is possible.
Although the system is capable of providing quantitative flow velocity and WSR measurement, further evaluation is required before clinical usage. A potential limitation exists in the idealized model, where only a straight and rigid tube is implemented. Therefore, a more complex geometry, such as the carotid bifurcation model previously simulated (Swillens et al. 2009), can be adapted. Such a framework, by combining the CFD with the ultrasound simulator, may provide a more realistic model. Further development should include vessel wall deformation, which has not been observed in current simulation. Another aspect for improvement of the simulation is the ultrasound image generation. The Field II ultrasound simulator is used in this study. It is a linear wave simulator that does not take into account the harmonic response of the simulated particles. This therefore restricts the modelling of the nonlinear response of microbubbles to generate contrast images similar to those in the in vitro experiment. However, the simulation of the contrast pulse sequence and response are possible and can be achieved by adaptation of the Creanius simulator (Varray et al. 2010).
Apart from the refinements in simulation, the limitations of using 2-D imaging to estimate velocity and WSR also need to be addressed. With 2-D imaging, only the inplane motion is tracked and considered when determining the accuracy of the developed system. However, the outof-plane motion may also affect the accuracy of wall velocity and shear rate measurement and therefore needs to be further investigated. For example, although a controlled flow rate is applied in the experiment, secondary flow and turbulence may arise because of the local tortuosity of the bifurcation model and may cause an underestimation in WSR. Despite this, the plane wave UIV remains capable of capturing the in-plane velocity and WSR distribution to show the relative difference between a healthy and a diseased carotid geometry.
In this work only the WSR has been calculated. The knowledge of fluid viscosity is required to calculate the WSS. Although blood is a non-Newtonian fluid, blood viscosity can be assumed as a constant in large arteries (Cheng et al. 2002;Kefayati et al. 2014). This is because the dependence of blood viscosity on shear rate is less pronounced in a large vessel than in the microcirculation. Such a phenomenon is known as the Fåhraeus-Lindqvist effect, which describes the variation of the hematocrit in relation to the vessel diameter (Fåhraeus and Lindqvist 1931). The WSS in small vessels, on the other hand, can also be calculated by taking into account the non-Newtonian characteristics of blood. The dynamic viscosity of the blood can either be measured directly using a viscometer (Jeong and Rosenson 2013) or fitted using a mathematical model (Soulis et al. 2008).
HFR plane wave UIV as a non-invasive technique is capable of evaluating both velocity and WSR simultaneously and potentially can be used in a wide range of clinical applications. For instance, the technique could help to predict atherosclerotic lesions, assess risk of thrombosis and provide insight into the pathogenesis to guide clinical decisions. It can also be used to predict the risk of angiographic restenosis after coronary and peripheral angioplasty and improve surgical techniques related to coronary artery bypass graft.

CONCLUSIONS
In this study, an improved plane wave UIV system to provide both flow velocity vectors and WSR mapping simultaneously has been developed and evaluated. Initial simulation studies have demonstrated the capabilities of this technique to provide accurate spatio-temporal varying velocity and WSR estimation. An in vitro experiment demonstrated the potential of this technique for vascular applications in vivo.