Advertisement

Spatio-Temporal Flow and Wall Shear Stress Mapping Based on Incoherent Ensemble-Correlation of Ultrafast Contrast Enhanced Ultrasound Images

Open AccessPublished:October 13, 2017DOI:https://doi.org/10.1016/j.ultrasmedbio.2017.08.930

      Abstract

      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.

      Key Words

      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 D.A.
      • Orekhov A.N.
      • Bobryshev Y.V.
      Effects of shear stress on endothelial cells: Go with the flow.
      ,
      • Chiu J.J.
      • Chien S.
      Effects of disturbed flow on vascular endothelium: Pathophysiological basis and clinical perspectives.
      ), but also trigger the biochemical event that results in the deposition and adherence of the platelet or the formation of atheroma (
      • Caro C.G.
      Discovery of the role of wall shear in atherosclerosis.
      ,
      • Cecchi E.
      • Giglioli C.
      • Valente S.
      • Lazzeri C.
      • Gensini G.F.
      • Abbate R.
      • Mannini L.
      Role of hemodynamic shear stress in cardiovascular disease.
      ,
      • Dhawan S.S.
      • Avati Nanjundappa R.P.
      • Branch J.R.
      • Taylor W.R.
      • Quyyumi A.A.
      • Jo H.
      • McDaniel M.C.
      • Suo J.
      • Giddens D.
      • Samady H.
      Shear stress and plaque development.
      ). 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 D.
      • Kaiktsis L.
      • Chaniotis A.
      • Pantos J.
      • Efstathopoulos E.P.
      • Marmarelis V.
      Wall shear stress: Theoretical considerations and methods of measurement.
      ,
      • Papaioannou T.G.
      • Stefanadis C.
      Vascular wall shear stress: Basic principles and methods.
      ). 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 J.P.
      • Wasserman B.A.
      • Steinman D.A.
      Errors in the estimation of wall shear stress by maximum Doppler velocity.
      ). 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 C.P.
      • Parker D.
      • Taylor C.A.
      Quantification of wall shear stress in large blood vessels using Lagrangian interpolation functions with cine phase-contrast magnetic resonance imaging.
      ,
      • Poelma C.
      • van der Mijle R.M.E.
      • Mari J.M.
      • Tang M.-X.
      • Weinberg P.D.
      • Westerweel J.
      Ultrasound imaging velocimetry: Toward reliable wall shear stress measurements.
      ). 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 A.M.
      • Frayne R.
      • Unal O.
      • Krupinski E.
      • Strother C.M.
      In vitro and in vivo comparison of three MR measurement methods for calculating vascular shear stress in the internal carotid artery.
      ,
      • Stalder A.F.
      • Russe M.F.
      • Frydrychowicz A.
      • Bock J.
      • Hennig J.
      • Markl M.
      Quantitative 2 D and 3 D phase contrast MRI: Optimized analysis of blood flow and vessel wall parameters.
      ) and ultrasound (
      • Poelma C.
      • van der Mijle R.M.E.
      • Mari J.M.
      • Tang M.-X.
      • Weinberg P.D.
      • Westerweel J.
      Ultrasound imaging velocimetry: Toward reliable wall shear stress measurements.
      ,
      • Zhang F.
      • Lanning C.
      • Mazzaro L.
      • Barker A.J.
      • Gates P.E.
      • Strain W.D.
      • Fulford J.
      • Gosling O.E.
      • Shore A.C.
      • Bellenger N.G.
      • Rech B.
      • Chen J.
      • Chen J.
      • Shandas R.
      In vitro and preliminary in vivo validation of echo particle image velocimetry in carotid vascular imaging.
      ). 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 C.H.
      • Bazigou E.
      • Eckersley R.J.
      • Yu A.C.H.
      • Weinberg P.D.
      • Tang M.X.
      Flow velocity mapping using high frame-rate plane wave ultrasound and microbubble contrast agents.
      ). 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:
      v(r)=v0(1r2R2)
      (1)


      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:
      τwall=μ2v0R
      (2)


      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 D.H.
      • McDicken W.N.
      Doppler ultrasound: Physics, instrumentation and signal processing.
      ). 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 J.R.
      Method for the calculation of velocity, rate of flow and viscous drag in arteries when the pressure gradient is known.
      ):
      Q=M10α2πR4μLΔPsin(ωt+ϕ+ε10)
      (3)


      where ΔPsin(ωt+ϕ) is the pressure difference, α is the Womersley number and M10 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 D.H.
      Some aspects of the relationship between instantaneous volumetric blood flow and continuous wave Doppler ultrasound recordings—III. The calculation of Doppler power spectra from mean velocity waveforms, and the results of processing these spectra with maximum, mean, and RMS frequency processors.
      ):
      vm(rR,t)=1πR2Qn|ψn|cos(ωtϕn+χn)
      (4)


      ψn(rR,τn)=τJ0(τn)τJ0(rRτn)τJ0(τn)2J1(τn)
      (5)


      where Jm is the mth order Bessel function of the first kind, τ=j3/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)=V0+n=1Vncos(nωtϕn)
      (6)


      The time evolution of the velocity profile can be determined by summing the velocity profile from each harmonic.
      V(rR,t)=2V0(1r2R2)+n=1Vn|ψn|cos(nωtϕn+χn)
      (7)


      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.
      Table 1Fourier component for flow velocities in the common carotid artery and common femoral artery
      Data adapted from Evans and McDicken (2000).
      Common carotid

      Diameter  = 6.0 mm

      Heart rate  = 60 bpm

      Viscosity  = 0.004 kgm–1 s–1
      Common femoral

      Diameter  = 8.4 mm

      Heart rate  = 60 bpm

      Viscosity  = 0.004 kgm–1 s–1
      nfαVnϕnnfαVnϕn
      01.0001.00
      11.033.90.337411.031.8932
      22.055.50.247922.057.72.4985
      33.086.80.2412133.089.51.28156
      44.107.80.1214644.1010.90.32193
      55.138.70.1114755.1312.20.27133
      66.159.60.1319766.1513.40.32155
      77.1810.30.0623377.1814.50.28195
      88.2112.40.0421888.2115.50.01310
      * Data adapted from
      • Evans D.H.
      • McDicken W.N.
      Doppler ultrasound: Physics, instrumentation and signal processing.
      .

      Ultrasound simulation

      Field II (
      • Jensen J.
      Field: A program for simulating ultrasound systems.
      ,
      • Jensen J.A.
      • Svendsen N.B.
      Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers.
      ) 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 P.R.
      The time-dependent force and radiation impedance on a piston in a rigid infinite planar baffle.
      ,
      • Tupholme G.E.
      Generation of acoustic pulses by baffled plane pistons.
      ). 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 J.M.
      Ultrasonic speckle formation, analysis and processing applied to tissue characterization.
      ). 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 J.A.
      • Svendsen N.B.
      Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers.
      ). In this case, the elements were further divided into four sub-elements 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 B.Y.
      • Tsang I.K.
      • Yu A.C.
      GPU-based beamformer: Fast realization of plane wave compounding and synthetic aperture imaging.
      ), in which delay and sum operations were applied to the RF data, was used to generate synthetic ultrasound images frame for further analysis.
      Table 2Field II simulation setup
      Probe parameterImaging parameter
      Centre frequency8 MHzImaging modePlane wave
      Number of element128Transmit frequency8 MHz
      Element pitch0.2 mmExcitation pulse1 cycle sinusoidal
      Element height5 mmPRF10 kHz
      Sampling frequency100 MHzCompounding angle5
      Elevational focus20 mmAngle Range20°
      Number of sub-aperture4
      PRF = pulse repetition frequency.

      In vitro flow experiment

      Ultrafast plane wave imaging system

      A research platform (Vantage 128, Verasonic,) was configured to acquire HFR plane wave ultrasound images. A linear 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 high-speed 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 C.H.
      • Bazigou E.
      • Eckersley R.J.
      • Yu A.C.H.
      • Weinberg P.D.
      • Tang M.X.
      Flow velocity mapping using high frame-rate plane wave ultrasound and microbubble contrast agents.
      ).

      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 P.S.
      • Luois S.
      • Dayton P.A.
      • Matsunaga T.O.
      Formulation and acoustic studies of a new phase-shift agent for diagnostic and therapeutic ultrasound.
      , resulting in a solution contained of 5 × 109 micobubbles/mL with an average size of 1 µm (
      • Sennoga C.A.
      • Yeh J.S.M.
      • Alter J.
      • Stride E.
      • Nihoyannopoulos P.
      • Seddon J.M.
      • Haskard D.O.
      • Hajnal J.V.
      • Tang M.-X.
      • Eckersley R.J.
      Evaluation of methods for sizing and counting of ultrasound contrast agents.
      ). 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 × 105 microbubbles/mL. This concentration is clinically relevant and has been used in previous studies (
      • Leow C.H.
      • Bazigou E.
      • Eckersley R.J.
      • Yu A.C.H.
      • Weinberg P.D.
      • Tang M.X.
      Flow velocity mapping using high frame-rate plane wave ultrasound and microbubble contrast agents.
      ,
      • Leow C.H.
      • Iori F.
      • Corbett R.
      • Duncan N.
      • Caro C.
      • Vincent P.
      • Tang M.X.
      Microbubble void imaging: A non-invasive technique for flow visualization and quantification of mixing in large vessels using plane wave ultrasound and controlled microbubble contrast agent destruction.
      ,
      • Tang M.-X.
      • Kamiyama N.
      • Eckersley R.J.
      Effects of nonlinear propagation in ultrasound contrast agent imaging.
      ).
      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 A.
      • Toulemonde M.
      • Yildiz Y.O.
      • Eckersley R.J.
      • Tang M.-X.
      Ultrasound imaging with microbubbles [Life Sciences].
      ). 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.
      Table 3Ultrafast data acquisition configurations
      Probe parameterImaging parameter
      ProbeL12-3 vImaging modePlane wave PI
      Centre frequency8 MHzTransmit frequency4 MHz
      Number of element128Excitation pulse1 cycle
      Element pitch0.2 mmPRF10 kHz
      Element height5 mmCompounding plane wave6 (3 positive, 3 negative)
      Sampling frequency32 MHzAngle range20° (10° step)
      Elevational focus20 mmImaging depth5 cm
      PI = pulse inversion; PRF = pulse repetition frequency.

      Velocity estimation using a modified-UIV analysis

      Advanced UIV analysis with an iterative window deformation model (
      • Scarano F.
      Iterative image deformation methods in PIV.
      ) was performed to estimate the flow velocity. Such a technique is implemented in a hierarchical scheme that progressively compensates the in-plane flow motion on a smaller scale. In brief, cross-correlations 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 J.
      • Scarano F.
      Universal outlier detection for PIV data.
      ) 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 S.T.
      • Meinhart C.D.
      Recent advances in micro-particle image velocimetry.
      ), also to increase the SNR of a correlation function to reduce correlation errors (
      • Meinhart C.D.
      • Wereley S.T.
      • Santiago J.G.
      A PIV algorithm for estimating time-averaged velocity fields.
      ,
      • Poelma C.
      • Mari J.M.
      • Foin N.
      • Tang M.-X.
      • Krams R.
      • Caro C.G.
      • Weinberg P.D.
      • Westerweel J.
      3D flow reconstruction using ultrasound PIV.
      ).
      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 B.
      • Tangen T.A.
      • Ekroll I.K.
      • Rolim N.
      • Torp H.
      • Bjåstad T.
      • Lovstakken L.
      Coherent plane wave compounding for very high frame rate ultrasonography of rapidly moving targets.
      ). 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.
      Fig. 1
      Fig. 1(a) Basic workflow of the modified UIV algorithm to generate velocity mapping overlaid on the HRIs. (b) Principle of incoherent ensemble-correlation approach where cross-correlation is applied to each a pair of LRIs to generate an averaged correlation with higher signal to noise ratio. UIV = ultrasound imaging velocity; HRIs = high-resolution images; LRIs = low-resolution images.

      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 A.
      • Golay M.J.
      Smoothing and differentiation of data by simplified least squares procedures.
      ). 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 S.
      • Tannenbaum A.
      localizing region-based active contours.
      ) 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:
      τw=με12
      (8)


      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 J.
      • Karri S.
      • Schmieg J.
      • Prabhu S.
      • Vlachos P.
      In vitro, time-resolved PIV comparison of the effect of stent design on wall shear stress.
      ,
      • Kundu P.K.
      Chapter 2—Cartesian tensors.
      ). The transformation can be written as
      εmn=i=12j=12CimCjnεij
      (9)


      where εij is the original strain tensor in the image Cartesian coordinate defined as
      εij=[uidxj+ujdxi]
      (10)


      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
      Cij=[cosθsinθsinθcosθ]
      (11)


      Fig. 2
      Fig. 2The 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.

      Error quantification

      The purpose of this study is to evaluate the performance of the developed technique, and two error metrics—the mean error (ME) and the normalized root mean square error (NRMSE)—were considered, as described in eqns (12) and (13).
      ME=0ty(t)dt0tyt(t)dt0yt(t)dt×100%
      (12)


      NRMSE=1yt(max)t=1n(yyt)2n×100%
      (13)


      where y and yt are the UIV estimated value, the reference value, yt(max) is the maximum reference value and n is the number of time sample.

      Anthropomorphic flow phantoms

      Two anatomically realistic flow phantoms (
      • Lai S.S.
      • Yiu B.Y.
      • Poon A.K.
      • Yu A.C.
      Design of anthropomorphic flow phantoms based on rapid prototyping of compliant vessel geometries.
      ), 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.

      Results

      Ultrasound flow simulation

      Steady flow

      The coherent and incoherent ensemble-correlation 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.
      Fig. 3
      Fig. 3Flow velocity vectors estimated using coherent (left) and incoherent ensemble correlation (right) under slow, medium and fast flow. The tube is tilted in an angle of 60° to the ultrasound beam.
      Table 4Velocity estimation accuracy of the coherent and incoherent ensemble-correlation approaches under different flow rates
      The error is calculated based on the radial velocity profile.
      V0 (cm/s)CoherentIncoherent
      ME (%)NRMSE (%)ME (%)NRMSE (%)
      20–6.08 ± 4.811.11 ± 0.03–6.01 ± 3.241.10 ± 0.04
      50–13.78 ± 3.753.35 ± 0.09–8.13 ± 2.452.68 ± 0.06
      80–69.16 ± 27.7412.96 ± 2.09–6.57 ± 1.545.55 ± 0.19
      ME = mean error; NRMSE = normalized root mean square error.
      * The error is calculated based on the radial velocity profile.
      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.
      Fig. 4
      Fig. 4(a) Vector velocity mapping derived from the simulated steady flow. (b–c) Comparison of the analytical reference value with the pre- and post-SG filter velocity profile and shear rate estimation. (d) Effects of filter length relative to the diameter of the tube and the filter order toward the WSR estimation error. SG = Savitsky-Golay; WSR = wall shear rate; n = filter length; dr = image resolution; D = vessel radial diameter.
      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 high-order polynomial is generally less robust and sensitive to small-amplitude and random distributed error (
      • Lou Z.
      • Yang W.J.
      • Stein P.D.
      Errors in the estimation of arterial wall shear rates that result from curve fitting of velocity profiles.
      ). 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 C.E.
      • Gharib M.
      Digital particle image velocimetry.
      ) and the non-slip boundary condition, which imposes a zero velocity at the boundary.
      Fig. 5
      Fig. 5(a) Effect of estimated wall location error toward the WSR estimation error. (b) Shear rate profile corresponding to various wall locations that deviate from the reference location. WSR = wall shear rate.

      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, Fig. 6, Fig. 7 show the temporal snapshot of the velocity vectors and WSR mapping at various phases of the pulsatile flow indicated in Fig. 6, Fig. 7. Spatial average 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.
      Fig. 6
      Fig. 6(a–b) Snapshots of the UIV estimation obtained at various phases in a simulated common carotid artery illustrate the key flow patterns and WSR mapping (left) and the comparison of the estimated flow profile with the analytical solution (right). (c) Comparison of the centerline velocity estimated using the UIV method and the analytical solution. The relative position of each flow pattern is marked in the centerline velocity plot and the shading in the UIV measurement represents the standard deviation. UIV = ultrasound imaging velocimetry; WSR = wall shear rate.
      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 velocity, while the direction of the near-wall flow velocity determines the direction in the WSR. A consistent WSR distribution across the field of view are illustrated; however, slight discrepancies can be observed near the left and right boundaries of the field of view or the edges in the lateral position as illustrated in Fig. 6, Fig. 7. This can be accounted for by the in-plane motion where the flow moves out of the field of view and the decreased spatial resolution in a compounded plane wave imaging.
      Fig. 7
      Fig. 7(a–e) Quantitative measurements of the UIV obtained at four phases in a simulated common femoral artery highlight the key flow patterns and WSR mapping (left) and the comparison of the estimated flow profile with the analytical solution (right). (e) Estimated centerline velocity overlaid on the analytical solution to show the accuracy of UIV. The relative position of each flow pattern is marked in the centerline velocity plot and the shading in the UIV measurement represents the standard deviation. UIV = ultrasound imaging velocimetry; WSR = wall shear rate.
      The centerline velocities and the spatial average velocity profiles are overlaid on the analytical solution, 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 spatial-averaged 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.
      Fig. 8
      Fig. 8Spatial-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.
      Table 5Error analysis on the estimated spatial average wall shear rate
      Error measuresCCACFA
      Upper wallLower wallUpper wallLower wall
      ME (%)6.54.97.210.6
      NRMSE (%)4.74.65.05.2
      CCA = common carotid artery; CFA = common femoral artery; ME = mean error; NRMSE = normalized root mean square error.

      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.
      Fig. 9
      Fig. 9Key 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.
      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 high-velocity 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 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 proximal 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 S.
      • Milner J.S.
      • Holdsworth D.W.
      • Poepping T.L.
      In vitro shear stress measurements using particle image velocimetry in a family of carotid artery models: Effect of stenosis severity, plaque eccentricity, and ulceration.
      ,
      • Yiu B.Y.
      • Lai S.S.
      • Yu A.C.
      Vector projectile imaging: Time-resolved dynamic visualization of complex flow patterns.
      ).
      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.
      Fig. 10
      Fig. 10The comparison of the flow velocity mapping at peak systole demonstrates the robustness of the incoherent ensemble correlation approach to track highly accelerated flow when a similar number of high-resolution images (N) are used to estimate the flow velocity. Note that the coherent ensemble correlations are estimated from the final compounded image (N) but the incoherent ensemble correlations are computed on the multiple low-resolution images (M × N) used to construct the same compounded image.

      Discussion

      The work presented in this study is an extension of previous work on an HFR-UIV system (
      • Leow C.H.
      • Bazigou E.
      • Eckersley R.J.
      • Yu A.C.H.
      • Weinberg P.D.
      • Tang M.X.
      Flow velocity mapping using high frame-rate plane wave ultrasound and microbubble contrast agents.
      ). 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 full-field of view velocity measurement is demonstrated in Fig. 4, Fig. 6, Fig. 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 B.
      • Fraser K.H.
      • Poelma C.
      • Mari J.-M.
      • Eckersley R.J.
      • Weinberg P.D.
      • Tang M.-X.
      Ultrasound imaging velocimetry: Effect of beam sweeping on velocity estimation.
      ).
      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 R.S.
      • Rittgers S.E.
      Derivation of shear rates from near-wall LDA measurements under steady and pulsatile flow conditions.
      ). 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 over-smoothing 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 Fig. 6, Fig. 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 over-estimation 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 H.-B.
      • Hertzberg J.
      • Lanning C.
      • Shandas R.
      Noninvasive measurement of steady and pulsating velocity profiles and shear rates in arteries using echo PIV: In vitro validation studies.
      ,
      • Masaryk A.M.
      • Frayne R.
      • Unal O.
      • Krupinski E.
      • Strother C.M.
      In vitro and in vivo comparison of three MR measurement methods for calculating vascular shear stress in the internal carotid artery.
      ). 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 contrast-specific 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=|20log10(SmicrobubbleSTissue)|, was calculated. Smicrobubble is the mean microbubbles signal and STissue 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.
      Fig. 11
      Fig. 11The 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.
      The capabilities of plane wave UIV to quantify and visualize complex flow dynamics within a physiologic-relevant 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 S.
      • Milner J.S.
      • Holdsworth D.W.
      • Poepping T.L.
      In vitro shear stress measurements using particle image velocimetry in a family of carotid artery models: Effect of stenosis severity, plaque eccentricity, and ulceration.
      ). 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-of-plane 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 Fig. 3, Fig. 10. 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 out-of-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 measurement 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 non-contrast HFR ultrasound imaging. The coherent ensemble-correlation approach, with an integrated motion correction technique to remove motion artefacts in the coherent compounded images, has been demonstrated in the speckle-tracking technique without contrast agents (
      • Joos P.
      • Porée J.
      • Liebgott H.
      • Vray D.
      • Cloutier G.
      • Nicolas B.
      • Garcia D.
      High-frame-rate velocity vector imaging echocardiography: An in vitro evaluation.
      ). 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 low- resolution 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 non-contrast 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 J.
      • Nikolov S.
      • Yu A.C.H.
      • Garcia D.
      Ultrasound vector flow imaging: II: Parallel systems.
      ,
      • Jensen J.A.
      • Nikolov S.
      • Yu A.C.H.
      • Garcia D.
      Ultrasound Vector Flow Imaging: I: Sequential systems.
      ). Although promising results to track two-dimensional 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 A.
      • Segers P.
      • Torp H.
      • Løvstakken L.
      Two-dimensional blood velocity estimation with ultrasound: Speckle tracking versus crossed-beam vector Doppler based on flow simulations in a carotid bifurcation model.
      ). Our approach, on the other hand, 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 real-time 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 A.
      • Løvstakken L.
      • Kips J.
      • Torp H.
      • Segers P.
      Ultrasound simulation of complex flow velocity fields based on computational fluid dynamics.
      ), 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 F.
      • Cachard C.
      • Tortoli P.
      • Basset O.
      Nonlinear radio frequency image simulation for harmonic imaging: Creanuis.
      ).
      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 in-plane motion is tracked and considered when determining the accuracy of the developed system. However, the out-of-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 C.P.
      • Parker D.
      • Taylor C.A.
      Quantification of wall shear stress in large blood vessels using Lagrangian interpolation functions with cine phase-contrast magnetic resonance imaging.
      ,
      • Kefayati S.
      • Milner J.S.
      • Holdsworth D.W.
      • Poepping T.L.
      In vitro shear stress measurements using particle image velocimetry in a family of carotid artery models: Effect of stenosis severity, plaque eccentricity, and ulceration.
      ). 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åhræus-Lindqvist effect, which describes the variation of the hematocrit in relation to the vessel diameter (
      • Fåhræus R.
      • Lindqvist T.
      The viscosity of the blood in narrow capillary tubes.
      ). 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 S.K.
      • Rosenson R.S.
      Shear rate specific blood viscosity and shear stress of carotid artery duplex ultrasonography in patients with lacunar infarction.
      ) or fitted using a mathematical model (
      • Soulis J.V.
      • Giannoglou G.D.
      • Chatzizisis Y.S.
      • Seralidou K.V.
      • Parcharidis G.E.
      • Louridas G.E.
      Non-Newtonian models for molecular viscosity and wall shear stress in a 3D reconstructed human left coronary artery.
      ).
      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.

      Acknowledgments

      Chee Hau Leow is supported by the Postgraduate Research Scholarship from Public Service Department of Malaysia . Meng-Xing Tang acknowledges the funding from the United Kingdom's Engineering and Physical Sciences Research Council (EPSRC) ( EP/M011933/1 ). The authors are grateful to Dr. Alfred Yu and his group for providing the carotid bifurcation phantoms and to Matthew Lee for the segmentation algorithm used in this report. The authors also thank Dr. Robert Eckersley and the late Professor David Cosgrove for the fruitful discussions.

      Supplementary Data

      The following is the supplementary data to this article:

      References

        • Caro C.G.
        Discovery of the role of wall shear in atherosclerosis.
        Arterioscler Thromb Vasc Biol. 2009; 29: 158-161
        • Cecchi E.
        • Giglioli C.
        • Valente S.
        • Lazzeri C.
        • Gensini G.F.
        • Abbate R.
        • Mannini L.
        Role of hemodynamic shear stress in cardiovascular disease.
        Atherosclerosis. 2011; 214: 249-256
        • Charonko J.
        • Karri S.
        • Schmieg J.
        • Prabhu S.
        • Vlachos P.
        In vitro, time-resolved PIV comparison of the effect of stent design on wall shear stress.
        Ann Biomed Eng. 2009; 37: 1310-1321
        • Cheng C.P.
        • Parker D.
        • Taylor C.A.
        Quantification of wall shear stress in large blood vessels using Lagrangian interpolation functions with cine phase-contrast magnetic resonance imaging.
        Ann Biomed Eng. 2002; 30: 1020-1032
        • Chistiakov D.A.
        • Orekhov A.N.
        • Bobryshev Y.V.
        Effects of shear stress on endothelial cells: Go with the flow.
        Acta Physiol (Oxf). 2017; 219: 382-408
        • Chiu J.J.
        • Chien S.
        Effects of disturbed flow on vascular endothelium: Pathophysiological basis and clinical perspectives.
        Physiol Rev. 2011; 91: 327-387
        • Denarie B.
        • Tangen T.A.
        • Ekroll I.K.
        • Rolim N.
        • Torp H.
        • Bjåstad T.
        • Lovstakken L.
        Coherent plane wave compounding for very high frame rate ultrasonography of rapidly moving targets.
        IEEE Trans Med Imaging. 2013; 32: 1265-1276
        • Dhawan S.S.
        • Avati Nanjundappa R.P.
        • Branch J.R.
        • Taylor W.R.
        • Quyyumi A.A.
        • Jo H.
        • McDaniel M.C.
        • Suo J.
        • Giddens D.
        • Samady H.
        Shear stress and plaque development.
        Expert Rev Cardiovasc Ther. 2010; 8: 545-556
        • Evans D.H.
        Some aspects of the relationship between instantaneous volumetric blood flow and continuous wave Doppler ultrasound recordings—III. The calculation of Doppler power spectra from mean velocity waveforms, and the results of processing these spectra with maximum, mean, and RMS frequency processors.
        Ultrasound Med Biol. 1982; 8: 617-623
        • Evans D.H.
        • McDicken W.N.
        Doppler ultrasound: Physics, instrumentation and signal processing.
        John Wiley & Sons, Hoboken, NJ2000
        • Fatemi R.S.
        • Rittgers S.E.
        Derivation of shear rates from near-wall LDA measurements under steady and pulsatile flow conditions.
        J Biomech Eng. 1994; 116: 361-368
        • Fåhræus R.
        • Lindqvist T.
        The viscosity of the blood in narrow capillary tubes.
        Am J Physiol Leg Content. 1931; 96: 562-568
        • Jensen J.
        Field: A program for simulating ultrasound systems.
        Med Biol Eng Comput. 1996; 34: 351-353
        • Jensen J.
        • Nikolov S.
        • Yu A.C.H.
        • Garcia D.
        Ultrasound vector flow imaging: II: Parallel systems.
        IEEE Trans Ultrason Ferroelectr Freq Control. 2016; 63: 1722-1732
        • Jensen J.A.
        • Nikolov S.
        • Yu A.C.H.
        • Garcia D.
        Ultrasound Vector Flow Imaging: I: Sequential systems.
        IEEE Trans Ultrason Ferroelectr Freq Control. 2016; 63: 1704-1721
        • Jensen J.A.
        • Svendsen N.B.
        Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers.
        IEEE Trans Ultrason Ferroelectr Freq Control. 1992; 39: 262-267
        • Jeong S.K.
        • Rosenson R.S.
        Shear rate specific blood viscosity and shear stress of carotid artery duplex ultrasonography in patients with lacunar infarction.
        BMC Neurol. 2013; 13: 36
        • Joos P.
        • Porée J.
        • Liebgott H.
        • Vray D.
        • Cloutier G.
        • Nicolas B.
        • Garcia D.
        High-frame-rate velocity vector imaging echocardiography: An in vitro evaluation.
        (In: 2016 IEEE International Ultrasonics Symposium (IUS). New York: IEEE, 1–4)2016
        • Katritsis D.
        • Kaiktsis L.
        • Chaniotis A.
        • Pantos J.
        • Efstathopoulos E.P.
        • Marmarelis V.
        Wall shear stress: Theoretical considerations and methods of measurement.
        Prog Cardiovasc Dis. 2007; 49: 307-329
        • Kefayati S.
        • Milner J.S.
        • Holdsworth D.W.
        • Poepping T.L.
        In vitro shear stress measurements using particle image velocimetry in a family of carotid artery models: Effect of stenosis severity, plaque eccentricity, and ulceration.
        PLoS ONE. 2014; 9 (e98209)
        • Kim H.-B.
        • Hertzberg J.
        • Lanning C.
        • Shandas R.
        Noninvasive measurement of steady and pulsating velocity profiles and shear rates in arteries using echo PIV: In vitro validation studies.
        Ann Biomed Eng. 2004; 32: 1067-1076
        • Kundu P.K.
        Chapter 2—Cartesian tensors.
        in: Cohen I.M. Dowling D.R. Fluid mechanics. 5th ed. Academic Press, Boston2012: 39-64
        • Lai S.S.
        • Yiu B.Y.
        • Poon A.K.
        • Yu A.C.
        Design of anthropomorphic flow phantoms based on rapid prototyping of compliant vessel geometries.
        Ultrasound Med Biol. 2013; 39: 1654-1664
        • Lankton S.
        • Tannenbaum A.
        localizing region-based active contours.
        IEEE Trans Image Process. 2008; 17: 2029-2039
        • Leow C.H.
        • Bazigou E.
        • Eckersley R.J.
        • Yu A.C.H.
        • Weinberg P.D.
        • Tang M.X.
        Flow velocity mapping using high frame-rate plane wave ultrasound and microbubble contrast agents.
        Ultrasound Med Biol. 2015; 41: 2913-2925
        • Leow C.H.
        • Iori F.
        • Corbett R.
        • Duncan N.
        • Caro C.
        • Vincent P.
        • Tang M.X.
        Microbubble void imaging: A non-invasive technique for flow visualization and quantification of mixing in large vessels using plane wave ultrasound and controlled microbubble contrast agent destruction.
        Ultrasound Med Biol. 2015; 41: 2926-2937
        • Lou Z.
        • Yang W.J.
        • Stein P.D.
        Errors in the estimation of arterial wall shear rates that result from curve fitting of velocity profiles.
        J Biomech. 1993; 26: 383-390
        • Masaryk A.M.
        • Frayne R.
        • Unal O.
        • Krupinski E.
        • Strother C.M.
        In vitro and in vivo comparison of three MR measurement methods for calculating vascular shear stress in the internal carotid artery.
        AJNR Am J Neuroradiol. 1999; 20: 237-245
        • Meinhart C.D.
        • Wereley S.T.
        • Santiago J.G.
        A PIV algorithm for estimating time-averaged velocity fields.
        J Fluids Eng. 2000; 122: 285-289
        • Mynard J.P.
        • Wasserman B.A.
        • Steinman D.A.
        Errors in the estimation of wall shear stress by maximum Doppler velocity.
        Atherosclerosis. 2013; 227: 259-266
        • Papaioannou T.G.
        • Stefanadis C.
        Vascular wall shear stress: Basic principles and methods.
        Hellenic J Cardiol. 2005; 46: 9-15
        • Poelma C.
        • Mari J.M.
        • Foin N.
        • Tang M.-X.
        • Krams R.
        • Caro C.G.
        • Weinberg P.D.
        • Westerweel J.
        3D flow reconstruction using ultrasound PIV.
        Exp Fluids. 2011; 50: 777-785
        • Poelma C.
        • van der Mijle R.M.E.
        • Mari J.M.
        • Tang M.-X.
        • Weinberg P.D.
        • Westerweel J.
        Ultrasound imaging velocimetry: Toward reliable wall shear stress measurements.
        Eur J Mech B Fluids. 2012; 35: 70-75
        • Savitzky A.
        • Golay M.J.
        Smoothing and differentiation of data by simplified least squares procedures.
        Anal Chem. 1964; 36: 1627-1639
        • Scarano F.
        Iterative image deformation methods in PIV.
        Meas Sci Technol. 2002; 13
        • Sennoga C.A.
        • Yeh J.S.M.
        • Alter J.
        • Stride E.
        • Nihoyannopoulos P.
        • Seddon J.M.
        • Haskard D.O.
        • Hajnal J.V.
        • Tang M.-X.
        • Eckersley R.J.
        Evaluation of methods for sizing and counting of ultrasound contrast agents.
        Ultrasound Med Biol. 2012; 38: 834-845
        • Sheeran P.S.
        • Luois S.
        • Dayton P.A.
        • Matsunaga T.O.
        Formulation and acoustic studies of a new phase-shift agent for diagnostic and therapeutic ultrasound.
        Langmuir. 2011; 27: 10412-10420
        • Soulis J.V.
        • Giannoglou G.D.
        • Chatzizisis Y.S.
        • Seralidou K.V.
        • Parcharidis G.E.
        • Louridas G.E.
        Non-Newtonian models for molecular viscosity and wall shear stress in a 3D reconstructed human left coronary artery.
        Med Eng Phys. 2008; 30: 9-19
        • Stalder A.F.
        • Russe M.F.
        • Frydrychowicz A.
        • Bock J.
        • Hennig J.
        • Markl M.
        Quantitative 2 D and 3 D phase contrast MRI: Optimized analysis of blood flow and vessel wall parameters.
        Magn Reson Med. 2008; 60: 1218-1231
        • Stanziola A.
        • Toulemonde M.
        • Yildiz Y.O.
        • Eckersley R.J.
        • Tang M.-X.
        Ultrasound imaging with microbubbles [Life Sciences].
        IEEE Signal Process Mag. 2016; 33: 111-117
        • Stepanishen P.R.
        The time-dependent force and radiation impedance on a piston in a rigid infinite planar baffle.
        J Acoust Soc Am. 1971; 49: 841-849
        • Swillens A.
        • Løvstakken L.
        • Kips J.
        • Torp H.
        • Segers P.
        Ultrasound simulation of complex flow velocity fields based on computational fluid dynamics.
        IEEE Trans Ultrason Ferroelectr Freq Control. 2009; 56: 546-556
        • Swillens A.
        • Segers P.
        • Torp H.
        • Løvstakken L.
        Two-dimensional blood velocity estimation with ultrasound: Speckle tracking versus crossed-beam vector Doppler based on flow simulations in a carotid bifurcation model.
        IEEE Trans Ultrason Ferroelectr Freq Control. 2010; 57: 327-339
        • Tang M.-X.
        • Kamiyama N.
        • Eckersley R.J.
        Effects of nonlinear propagation in ultrasound contrast agent imaging.
        Ultrasound Med Biol. 2010; 36: 459-466
        • Thijssen J.M.
        Ultrasonic speckle formation, analysis and processing applied to tissue characterization.
        Pattern Recognit Lett. 2003; 24: 659-675
        • Tupholme G.E.
        Generation of acoustic pulses by baffled plane pistons.
        Mathematika. 1969; 16: 209-224
        • Varray F.
        • Cachard C.
        • Tortoli P.
        • Basset O.
        Nonlinear radio frequency image simulation for harmonic imaging: Creanuis.
        (In: IEEE Ultrasonics Symposium. New York: IEEE, 2179–2182)2010
        • Wereley S.T.
        • Meinhart C.D.
        Recent advances in micro-particle image velocimetry.
        Annu Rev Fluid Mech. 2010; 42: 557-576
        • Westerweel J.
        • Scarano F.
        Universal outlier detection for PIV data.
        Exp Fluids. 2005; 39: 1096-1100
        • Willert C.E.
        • Gharib M.
        Digital particle image velocimetry.
        Exp Fluids. 1991; 10: 181-193
        • Womersley J.R.
        Method for the calculation of velocity, rate of flow and viscous drag in arteries when the pressure gradient is known.
        J Physiol. 1955; 127: 553-563
        • Yiu B.Y.
        • Lai S.S.
        • Yu A.C.
        Vector projectile imaging: Time-resolved dynamic visualization of complex flow patterns.
        Ultrasound Med Biol. 2014; 40: 2295-2309
        • Yiu B.Y.
        • Tsang I.K.
        • Yu A.C.
        GPU-based beamformer: Fast realization of plane wave compounding and synthetic aperture imaging.
        IEEE Trans Ultrason Ferroelectr Freq Control. 2011; 58: 1698-1705
        • Zhang F.
        • Lanning C.
        • Mazzaro L.
        • Barker A.J.
        • Gates P.E.
        • Strain W.D.
        • Fulford J.
        • Gosling O.E.
        • Shore A.C.
        • Bellenger N.G.
        • Rech B.
        • Chen J.
        • Chen J.
        • Shandas R.
        In vitro and preliminary in vivo validation of echo particle image velocimetry in carotid vascular imaging.
        Ultrasound Med Biol. 2011; 37: 450-464
        • Zhou B.
        • Fraser K.H.
        • Poelma C.
        • Mari J.-M.
        • Eckersley R.J.
        • Weinberg P.D.
        • Tang M.-X.
        Ultrasound imaging velocimetry: Effect of beam sweeping on velocity estimation.
        Ultrasound Med Biol. 2013; 39: 1672-1681