Robust flow vector estimation for echocardiography with extended Nyquist velocity using dual-PRF approach: a flow phantom study

The Nyquist velocity ${v}_{N}$ in the color Doppler method is given by Eq. (1):

$$\begin{eqnarray}&&\begin{array}{c}{v}_{N}=\displaystyle \frac{{\rm{P}}{\rm{R}}{\rm{F}}\times \lambda }{4},\end{array}\end{eqnarray} \tag{ 1 }$$

where $\lambda$ is the wavelength of the acoustic pluses. In the dual-PRF method, two different PRFs (PRF₁ and PRF₂) are selected, where PRF₁ is defined based on the imaging depth and expected maximum flow velocity, while PRF₂ is set as shown in Eq. (2):

$$\begin{eqnarray}&&\begin{array}{c}{\rm{P}}{\rm{R}}{{\rm{F}}}_{2}=\displaystyle \frac{p}{q}{\rm{P}}{\rm{R}}{{\rm{F}}}_{1},\end{array}\end{eqnarray} \tag{ 2 }$$

where $p$ and $q$ are mutually prime numbers and $p\lt q.$ In addition, the relationship between the Nyquist speeds, ${v}_{{N}_{1}}$ and ${v}_{{N}_{2}}$, can be derived as follows:

$$\begin{eqnarray}&&\begin{array}{c}{v}_{{N}_{2}}=\displaystyle \frac{p}{q}{v}_{{N}_{1}}.\end{array}\end{eqnarray} \tag{ 3 }$$

In this study, we used a condition, $p=2$ and $q=3.$ As illustrated in Fig. 1, the Doppler shifts measured with the respective PRFs aliased at their respective Nyquist speeds and showed the relationships between the actual and measured velocities. Each color-coded region in Fig. 1 shows a unique difference between the flow speeds measured by ${\rm{P}}{\rm{R}}{{\rm{F}}}_{1}$ and ${\rm{P}}{\rm{R}}{{\rm{F}}}_{2}$ and can be identified using the Eq. (4):

$$\begin{eqnarray}&&\begin{array}{l}{D}_{v}=\displaystyle \frac{3\left({v}_{{D}_{2}}-{v}_{{D}_{1}}\right)}{2{v}_{{N}_{1}}},\\ {C}_{{\rm{a}}{\rm{l}}{\rm{i}}{\rm{a}}{\rm{s}}}={\rm{n}}{\rm{i}}{\rm{n}}{\rm{t}}\left({D}_{v}\right),\end{array}\end{eqnarray} \tag{ 4 }$$

where ${v}_{{D}_{1}}\,$and ${v}_{{D}_{2}}\,$are the flow speeds measured at the respective PRFs, and the ${\rm{n}}{\rm{i}}{\rm{n}}{\rm{t}}$ function gives the closest integer value. The value of ${C}_{{\rm{a}}{\rm{l}}{\rm{i}}{\rm{a}}{\rm{s}}}$ calculated in Eq. (4) was used to obtain the number of aliased cycles (Nyquist numbers, ${n}_{{N}_{i}}$) at each pixel in the color Doppler frames measured with the respective PRFs, as listed in Table I. Finally, the Nyquist numbers, ${n}_{{N}_{1}\,}\,$and ${n}_{{N}_{2}}$, were substituted into Eq. (5) to calculate the de-aliased Doppler velocity, $\overline{{v}_{D}}:$

$$\begin{eqnarray}&&\begin{array}{c}\overline{{v}_{D}}=\displaystyle \frac{2\left({v}_{{D}_{1}}+2{n}_{{N}_{1}}{v}_{{N}_{1}}\right)+3\left({v}_{{D}_{2}}+2{n}_{{N}_{2}}{v}_{{N}_{2}}\right)}{5}.\end{array}\end{eqnarray} \tag{ 5 }$$

Zoom In Zoom Out Reset image size

Table I. The Nyquist numbers (number of aliasing cycles for each PRF measurement) determined by different ${C}_{{\rm{a}}{\rm{l}}{\rm{i}}{\rm{a}}{\rm{s}}}$ values.

${{\boldsymbol{C}}}_{{\bf{alias}}}$	${{\boldsymbol{n}}}_{1}$	${{\boldsymbol{n}}}_{2}$
2	0	−1
1	1	1
0	0	0
−1	−1	−1
−2	0	1

As shown in Fig. 1, the Nyquist speed of the dual-PRF color Doppler could be increased to twice the original Nyquist speed, ${v}_{{N}_{1}}.$

After the dual-PRF de-aliasing was performed, aliasing-free color Doppler frames for each transmission angle were obtained. At each pixel in those frames, the flow speed components along each transmission angle, ${v}_{1}$ and ${v}_{2},$ were used to compute the 2D flow vectors. Fig. 2(a) illustrates the geometric relationships between the transmitted diverging waves and the measured and de-aliased flow speed components. The actual 2D velocity components along the x and z axes, ${V}_{x}$ and ${V}_{z},$ could be decomposed as shown in Figs. 2(b) and 2(c) using ${v}_{1}$ and ${v}_{2},$ and those relationships can be expressed as Eqs. (6) and (7):

$$\begin{eqnarray}&&\begin{array}{c}-{v}_{1}={V}_{x}\,\sin \,{\theta }_{1}+{V}_{z}\,\cos \,{\theta }_{1},\end{array}\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\begin{array}{c}-{v}_{2}={V}_{x}\,\sin \,{\theta }_{2}+{V}_{z}\,\cos \,{\theta }_{2},\end{array}\end{eqnarray} \tag{ 7 }$$

where ${\theta }_{1}$ and ${\theta }_{2}$ are the transmission angles and should be calculated at each pixel because of the nature of the diverging wave. Finally, from Eqs. (6) and (7), the components of the 2D velocity vector ($\widehat{{V}_{x}},$ $\widehat{{V}_{z}}$) are derived as follows:

$$\begin{eqnarray}&&\begin{array}{c}\widehat{{V}_{x}}=-\displaystyle \frac{{v}_{1}\,\cos \,{\theta }_{2}-{v}_{2}\,\cos \,{\theta }_{1}}{\sin \left({\theta }_{1}-{\theta }_{2}\right)},\end{array}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\begin{array}{c}\widehat{{V}_{z}}=\displaystyle \frac{{v}_{1}\,\sin \,{\theta }_{2}-{v}_{2}\,\sin \,{\theta }_{1}}{\sin \left({\theta }_{1}-{\theta }_{2}\right)}.\end{array}\end{eqnarray} \tag{ 9 }$$

Zoom In Zoom Out Reset image size

First, we evaluated the de-aliasing performance of the proposed dual-PRF, dual-angle flow measurement scheme because, as discussed in Sect. 2.2, it was important to accurately estimate the 1D flow speed in the color Doppler method to achieve accurate and robust estimations of the 2D velocities. The imaging system was comprised of a research-purpose ultrasound system (Vantage 256 system, Verasonics Inc., USA) equipped with a custom sector probe (64 channels, center frequency: 2.5 MHz). The flow phantom (Doppler 403 Flow Phantom, Model 1425B, GAMMEX, USA) had a straight flow channel with an inner diameter of 5 mm sloped at $40^\circ$ (Fig. 3). Two flow conditions were used (a constant flow of 8.0 ml s⁻¹ and pulsatile flow of 12 ml s⁻¹ at 60 bpm). In this study, the flow signals obtained from the flow channel at a depth of 5–8 cm (indicated as the red region in Fig. 3) were processed and visualized by manually setting the region of interest (RoI) based on a power Doppler image.

Zoom In Zoom Out Reset image size

We designed an imaging sequence for the dual-PRF, dual-angle scheme (Fig. 4) using the detailed imaging parameters listed in Table II. First, a B-mode dataset was acquired by emitting 15 diverging wave pulses (Tx angle range: $-42^\circ$–$+42^\circ$). Second, the dual-PRF Doppler sequence was performed as indicated by the yellow and green boxes in Fig. 4. In each block, thirty-four diverging wave pulses were alternatively emitted from $\pm 30^\circ .$ The PRF for the first and third blocks was set at 2.0 kHz $\left({\rm{P}}{\rm{R}}{{\rm{F}}}_{1}\right),$ while that of the other blocks was at 1.3 kHz $\left({\rm{P}}{\rm{R}}{{\rm{F}}}_{2}\right).$ In each data acquisition, 10 frames of data (B-mode and flow vector) were acquired for the constant flow measurements, and 100 frames were acquired for the pulsatile flow measurements.

Zoom In Zoom Out Reset image size

Table II. Ultrasound imaging parameters for the dual-PRF, dual-angle imaging framework.

Parameters	Values
Probe	Sector
Center frequency [MHz]	2.5
Channel number	64
B-mode
PRF [kHz]	4.5
Number of transmissions	15
Range of transmission angle [°]	$-42\mbox{--}+42$
Pulse length	1 pulse
Doppler
PRF1 [kHz]	2.0
Nyquist speed ${v}_{{N}_{1}}$ [mm s−1]	310
PRF2 [kHz]	1.3
Nyquist speed ${v}_{{N}_{2}}$ [mm s−1]	207
Extended Nyquist speed [mm s−1]	620
Pulse length	3 pluses
Ensemble size	65
Transmission angles [°]	$\pm 30$
Frame rate [fps]	20
Imaging depth [cm]	13.5

The acquired data were processed in a pipeline, as shown in Fig. 5. For the Doppler measurements at each PRF setting, the acquired echo data were processed using a singular value decomposition clutter filter and a lag-1 autocorrelator to generate color Doppler images. ^25–27) The aliasing artifacts in the color Doppler images obtained with the different transmission angles were corrected with the dual-PRF de-aliasing framework described in Sect. 2.1 using data obtained with two PRF conditions. To compare the de-aliasing performance of the dual-PRF method with that of the conventional manual phase unwrapping approach (hereafter, the single-PRF method), we used color Doppler data acquired at ${\rm{P}}{\rm{R}}{{\rm{F}}}_{1}$ and manually de-aliased by applying Eq. (10) to all of the negative speed values:

$$\begin{eqnarray}&&\begin{array}{c}\overline{{v}_{{D}_{1}}}={v}_{{D}_{1}}+2{v}_{{N}_{1}}.\end{array}\end{eqnarray} \tag{ 10 }$$

Zoom In Zoom Out Reset image size

After performing the respective de-aliasing methods, the estimated flow speed maps were processed with a temporal (9 frames) and spatial (6 $\times$ 6 pixels) smoothing filter.

To verify the accuracy of the 2D velocity vector estimation using the proposed imaging framework, another flow phantom experiment was conducted to measure the same flow with ultrasound and a PIV system (Kato Koken Co., Ltd., Japan). Figure 6 illustrates the setup used in this experiment. The straight flow channel phantom had an inner diameter of 8 mm and was made of PVA-gel. ^28–30) It is important to note that this PVA-gel phantom was both acoustically and optically transparent. Thus, the ultrasound and PIV methods could measure the same flow profile. The blood mimicking fluid (BMF) used as the liquid in the flow phantom consisted of distilled water, spherical silica particles (Godd ball, Suzukiyushi Industrial Co., Ltd., Japan) as an acoustic scatterer, and polyamide seeding particles (PSP-50, DANTEC, Denmark) as a PIV tracer. To circulate the BMF, the phantom was placed in a water tank and connected to a custom gear-flow pump to generate a constant flow of 18 ml s⁻¹ or a pulsatile flow of 21 ml s⁻¹ at 60 bpm. For ultrasound imaging, a sector probe was placed above the phantom at a tilt angle of $60^\circ$ [Fig. 6(a)]. For the PIV measurement, a laser sheet and high-speed camera were used to acquire the flow profile on the same imaging plane as the ultrasound imaging [Fig. 6(b)]. The high-speed camera captured the motion of the tracer particles at 1000 fps, and the 2D flow vector field was calculated using commercial PIV software (Flow Expert 2D2C, Kato Koken).

Zoom In Zoom Out Reset image size

In this experiment, we used PRF values (3.0 kHz for PRF₁ and 2.0 kHz for PRF₂) to acquire images at a depth of approx. 10.2 cm (Table III) and processed the acquired echo data with the data processing pipeline shown in Fig. 5. In addition to the dual-angle, dual-PRF flow velocity estimation, we also performed 2D flow velocity estimation using the dual-angle color Doppler data measured at PRF₁ without any aliasing correction to verify the efficacy of the proposed method.

Table III. PRFs and their Nyquist speed values used for the experiments to compare the 2D flow vectors estimated with the proposed and PIV measurements.

Parameters	Values
PRF1 [kHz]	3.0
Nyquist speed ${{\boldsymbol{v}}}_{{{\boldsymbol{N}}}_{1}}$ [mm s−1]	465
PRF2 [kHz]	2.0
Nyquist speed ${{\boldsymbol{v}}}_{{{\boldsymbol{N}}}_{2}}$ [mm s−1]	310
Extended Nyquist speed [mm s−1]	930
Pulse length	3 pluses
Ensemble size	65
Transmission angle [°]	$\pm 30$
Frame rate [fps]	27
Imaging depth [cm]	10.2

Figures 7(a) and 7(b) show examples of the obtained color Doppler images (transmission angle: −30°) under the different PRF conditions in the constant flow setting. Because of the variation in the diverging beam direction [Fig. 7(c)], the measured flow speed and appearance of the aliasing artifacts in the color Doppler images changed with the depth. For example, the beam directions at the red and white points in Fig. 7(c) were 2.4° and 17.5°, respectively, and the angles between the beam direction and the flow direction (40°) were 47.6° and 32.5°, respectively. As such, in this imaging scenario, the measured flow speed continuously increased from the upper part to the bottom part in the flow region, and aliasing artifacts occurred when the measured flow speed exceeded the Nyquist velocity of respective PRFs [Figs. 7(a) and 7(b)]. Using those color Doppler images, aliasing-free color Doppler images were generated based on the dual-PRF de-aliasing approach [Fig. 7(d)] and showed profiles similar to those obtained with the single-PRF method [Fig. 7(e)]. The percentage errors between the images de-aliased with the dual-PRF and single-PRF methods [Fig. 7(f)] also showed high similarity between the two methods, and the mean absolute percentage error was found to be 21.5%. In the upper right region of the flow tract, where the angle between the acoustic beam and flow direction was high and the measured flow speed was small, it was observed that the errors between the dual-PRF and single-PRF methods were quite small, whereas in the lower left region of the tract, where the measured flow speed was high, the error values were larger.

Zoom In Zoom Out Reset image size

The de-aliasing performance of the proposed framework was also evaluated in the pulsatile flow condition. Figures 8(a)–8(d) shows the color Doppler images (transmission angle of $-30^\circ$) acquired at slow and fast flow phases in a cycle of pulsatile flow. Similar to the constant flow images, a higher flow speed was measured in the lower left region, corresponding to the beam direction [Fig. 8(k)]. In the fast flow phase, color Doppler images in both PRF conditions [Figs. 8(b) and 8(d)] had aliasing artifacts. Particularly in PRF₂ condition [Fig. 8(d)], the flow velocity was measured to be positive without aliasing in the upper right area, aliased to the negative with one-cycle of aliasing in the middle area, and finally increased to the positive value with one-cycle of aliasing in the lower left area. This observation could also be found in the flow velocity profile shown in Fig. 8(l). With this setting, the dual-PRF method effectively corrected the aliasing artifacts in the fast flow phase [Fig. 8(f)]. In addition, it maintained the appearance of the flow profile observed in the slow flow phase [Fig. 8(e)], where the original color Doppler image [Figs. 8(a) and 8(c)] did not have any aliasing artifacts. Figures 8(i) and 8(j) show the respective percentage errors between the images obtained with the dual-PRF and single-PRF methods. The mean percent error was 3.9% in the slow flow phase and increased to 32.5% in the fast flow phase. In the fast flow phase, large errors were observed near the wall of the flow tract in the lower left area.

Zoom In Zoom Out Reset image size

3.1.1. Causes of errors observed with the dual-PRF de-aliasing method

Both the constant and pulsatile flow experiments showed that the dual-PRF approach successfully increased the Nyquist speed of the color Doppler to twice that under the PRF₁ condition, and the aliased regions colored in blue in Figs. 7(a) and 8(b) were appropriately corrected as shown in Figs. 7(d) and 8(f), respectively. As shown in Figs. 7(f) and 8(j), the de-aliased flow speed values were comparable to those calculated with the single-PRF approach (manual de-aliasing). However, high errors were found near the wall under the lower left area in both flow conditions, except for the slow flow phase of the pulsatile setting. It is suggested that the error became higher near the wall when the angle between the beam and flow directions became small. The following two factors might be related to those observations: (1) inadequate clutter filtering and (2) high variance in the Doppler signal at the boundary of the aliasing regions. When the clutter filtering was not appropriate (e.g. the filter excessively suppressed the slow flow velocities near the wall in the fast flow phase), the flow speed estimation in the color Doppler was affected and might have been underestimated, as shown in the expanded view of Fig. 7(b). The optimal clutter filtering approach for the proposed framework needs to be investigated in future study, e.g. using an adaptive adjustment of filtering parameters. ³¹⁾ Moreover, the gradient of the flow speed in the flow tract was observed to increase in the lower left area as the angle between the acoustic beam and low direction decreased. When the original color Doppler images had aliasing artifacts, the spatio-temporal variance in the Doppler signals at those regions also increased and might affect the flow speed estimation. ³²⁾ For instance, the color Doppler image measured at PRF₂ in the fast flow phase [Fig. 8 (d)] had a region of positive values with one-cycle of aliasing in the lower left area and might have a high-velocity gradient between the wall side and center of the flow tract.

In the dual-PRF de-aliasing approach, the errors in the original color Doppler images could affect the accuracy of the flow speed estimation because the method utilized the differences between the flow speed values measured at different PRFs. As shown in Eq. (4), the difference between two values (${D}_{v}$) was rounded using the ${\rm{n}}{\rm{i}}{\rm{n}}{\rm{t}}$ function and used as ${C}_{{\rm{a}}{\rm{l}}{\rm{i}}{\rm{a}}{\rm{s}}}$ to determine the aliasing conditions at each pixel. For example, under the settings for this experiment (${v}_{{N}_{1}}$ = 310 mm s⁻¹ at PRF₁: 2 kHz), the absolute difference between ${v}_{{D}_{1}}\,$and ${v}_{{D}_{2}}$ had to be in $\pm 100$ mm s⁻¹ of the expected difference shown in Fig. 1 or it would provide an incorrect ${C}_{{\rm{a}}{\rm{l}}{\rm{i}}{\rm{a}}{\rm{s}}}$ value. Figure 9(a) shows a map of the ${C}_{{\rm{a}}{\rm{l}}{\rm{i}}{\rm{a}}{\rm{s}}}$ values for the images in Figs. 8(b) and 8(d). The ${C}_{{\rm{a}}{\rm{l}}{\rm{i}}{\rm{a}}{\rm{s}}}$ value was properly estimated as 0 in the upper right area, which sequentially turns to −2 and then to 1 along the flow channel center, which corresponded to the observations in Fig. 8(l). However, the estimation of ${C}_{{\rm{a}}{\rm{l}}{\rm{i}}{\rm{a}}{\rm{s}}}$ value tended to be uncertain in the lower left region and the near-wall region, where $\left|{D}_{v}-{C}_{{\rm{a}}{\rm{l}}{\rm{i}}{\rm{a}}{\rm{s}}}\right|$ values were higher than 0.25 [Fig. 9(b)]. This uncertain region agreed with the high-error region, where the relative percent error in Fig. 8(j) is lower than −10% [Fig. 9(c)]. Thus, the dual-PRF approach might produce a spurious flow speed estimate when the original color Doppler has erroneous values. To improve the accuracy and robustness of aliasing-free flow speed estimation in the dual-PRF approach, it might be effective to combine the proposed method with a block matching scheme, ³³⁾ referring to the extended least-squares method previously reported. ³⁴⁾

Zoom In Zoom Out Reset image size

Figure 10(a) shows the flow vectors under the constant flow setting visualized using the dual-PRF, dual-angle vector Doppler method. It clearly shows that the method achieved a robust flow vector estimation that successfully corrected the aliasing artifacts, which could be seen as the wrong direction and speed [Fig. 10(b)] when the flow vectors were calculated without any aliasing correction. The flow vectors obtained using the proposed method were comparable to those measured with PIV [the ground-truth in this study; Fig. 10(c)]. In addition, the flow speed profiles measured with the proposed method and PIV are shown in Fig. 11. For each profile, the absolute value of the estimated flow vectors along the short-axis cross section of the flow tract was calculated and plotted. The profile data were generated by averaging the 2D flow vectors over 100 frames and approximately 11 mm along the flow direction from the center of the image. Although the dual-PRF de-aliasing method showed a spatial flow speed profile that was similar to that of the PIV method, the peak flow velocity obtained with dual-PRF, dual-angle vector Doppler was lower (627.4 $\pm \,$42.7 mm s⁻¹) than the value of 849.3 $\pm \,$3.9 mm s⁻¹ measured by the PIV method, and the relative error at the flow channel center was −26.1%. The error might have occurred because the flow vector estimations based on the color Doppler imaging and least-squares fitting methods tended to underestimate the peak flow speed as a result of the nature of the mean phase-shift estimations. ^19,35)

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Furthermore, the flow vectors under the pulsatile flow setting (for both the slow and fast flow phases) measured with the proposed method and PIV are shown in Fig. 12, and their flow speed profiles along the short-axis cross section are shown in Fig. 13. The dual-PRF and dual-angle vector Doppler methods yielded flow velocity profiles similar to those of the PIV method. As also shown in Fig. 11, the estimated peak flow speed (605.9 $\pm \,$106.5 mm s⁻¹) obtained with the proposed method was slightly lower than that (708.5 $\pm \,$8.6 mm s⁻¹) obtained with the PIV method in the fast flow phase, while the proposed method showed a slightly higher speed (162.1 $\pm \,$24.6 mm s⁻¹) than that of PIV (119.9 $\pm \,$5.2 mm s⁻¹) in the slow flow phase. The relative percent error at the flow channel center were −14.5% in the fast flow phase and 35.2% in the slow flow phase, respectively. This observation might be related to the spatial and temporal variances of the flow speed, which would have affected the color Doppler measurements. ³⁵⁾

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Both in the constant and pulsatile conditions, the vector estimation was observed to be unstable especially at the lower left area and near-wall region [Figs. 10(a), 12(b)]. There are two possible reasons on this point: (1) higher color Doppler errors occurred in the de-aliasing process as mentioned in Sect. 3.1.1, and (2) decreased differences between two transmission angles. Due to the nature of dual-angle diverging wave transmission, the difference between the two transmission angles becomes smaller in the deeper region that affects the accuracy of the vector estimation as shown in the previous studies. ¹⁹⁾ Moreover, either in the proposed method or PIV, the spatial flow speed profiles (Figs. 11 and 13) were not observed to be parabolic as those expected in the straight tube phantom. This might be caused by multiple factors such as the deformation and surface properties of the flow tract and turbulent flow component as the Reynolds number was approximately 3000. ³⁶⁾ Further evaluation of the accuracy of the proposed method is necessary, such as by precisely synchronizing two measurement modalities that will enable simultaneous acquisition of the target flow profiles.

Finally, with the imaging parameters used in this study, the theoretical extended Nyquist speed was approx. 930 mm s⁻¹ as listed in Table III, which was slightly lower than the peak flow speed (approx. 1000 mm s⁻¹) in the left ventricle. This study only tested the proposed framework with straight tube models, and the maximum flow speed did not reach the extended Nyquist flow speed limit. Also, the present study has not taken into account complex flow patterns such as vortex flow and flow jet that would present in the heart. Therefore, in future studies, the efficacy of the proposed method and the optimal imaging parameters including the combination of p and q in Eq. (2) and the PRFs, should be investigated further in more realistic models, such as a left ventricle model. ³⁷⁾