A Color-Doppler Shear-Wave-Imaging Phase-reconstruction Method Using Four Color Flow Images

Abstract

This study investigates shear wave phase map reconstruction using a limited number of color flow images (CFIs) acquired with a color Doppler ultrasound imaging instrument. We propose an efficient reconstruction method to considerably reduce the number of CFIs required for reconstruction and compare this method with Fourier analysis-based color Doppler shear wave imaging. The proposed method uses a two-step phase reconstruction process, including an initial phase map derived from four CFIs using an advanced iterative algorithm of optical interferometry. The second step reduces phase artifacts in the initial phase map using an iterative correction procedure that cycles between the Fourier and inverse Fourier domains while imposing directional filtering and total variation regularization. We demonstrate the efficacy of this method using synthetic and experimental data of a breast phantom and human breast tissue. Our results show that the proposed method maintains image quality and reduces the number of CFIs required to four; previous methods have required at least 32 CFIs to achieve equivalent image quality. The proposed method is applicable to real-time shear wave elastography using a continuous shear wave produced by a mechanical vibrator.

Keywords

color Doppler shear wave imaging phase reconstruction iterative method total variation shear wave elasticity imaging

Introduction

Shear wave elasticity imaging (SWEI) delineates the elastic properties of soft tissues. It is used clinically to detect tumors and determine the characteristics and stiffness of patients’ organs.^1-4 Several methods for SWEI have been proposed, including strain elastography and ultrasound push pulse-based shear wave elastography.^5-11 Strain elastography has a number of advantages over other elastography methods, including quick and easy image acquisition. However, because quantification using this method is difficult, relative evaluations, such as the Tsukuba-score for thyroid nodules, are generally employed. Conversely, shear wave elastography enables quantitative information regarding soft tissue target areas to be acquired only if the shear wave propagates into the tissue. Ultrasound push pulse-based shear wave elastography uses acoustic radiation force impulses (ARFI) to generate shear wave propagations. Continuous or frequent excitations cannot be used because ARFI excitation can damage tissue through heat generation.

In 2015, Yamakoshi et al. developed a novel color Doppler shear wave imaging (CD SWI) technique.^10,11 In CD SWI, shear wave propagations are continuously generated using a mechanical vibrator on the body’s surface. Furthermore, CD SWI uses the color Doppler imaging method to acquire the flow information in vessels without contrast agents. Color Doppler imaging is one of the most commonly employed clinical ultrasound imaging functions; consequently, CD SWI can be performed using many routine ultrasound imaging systems.

In CD SWI, the wavefront of the shear wave is displayed as color flow images (CFIs) when the shear wave frequency and displacement amplitude conditions are met.¹⁰ The shear wave frequency condition is related to the pulse repetition frequency (PRF) of the ultrasonic imaging system, and the shear wave vibration frequency is represented by f_b = (k + 1 / 2) / 2Δt, where k is zero or an integer value and Δt is the time duration between successive ultrasonic pulses that radiate in the same direction for color flow velocity estimation, that is, the reciprocal of PRF. The shear wave displacement amplitude condition is related to the wavelength of the ultrasonic wave, which defines the minimum shear wave displacement amplitude required for image reconstruction, and the shear wave displacement amplitude is represented by 1/8 λ < ξ₀ < 3/8 λ, where λ is the wavelength of the ultrasonic wave. For example, when λ, k, and PRF are 230 µm, 1, and 365 Hz, 1/8 λ and f_b are approximately 30 µm and 274 Hz, which can easily be produced by a small mechanical vibrator.

The propagation of the wavefront of the shear wave, that is, its motion, can be visualized on the CFI by slightly altering the shear wave frequency from the frequency condition. Because the wave motion has a periodic function in the time domain, the shear wave phase can be acquired as a two-dimensional (2D) map using frequency analysis. Shear wave velocity and phase propagation direction maps can be generated from the shear wave phase maps.

CFIs including the shear wave wavefront can contain various artifacts from standing waves or reflected waves at tissue boundaries. To suppress the phase variations on shear wave phase maps, the current Fourier analysis requires at least 32 CFI frames to be captured every 2 to 3 s.¹¹ The phase amount on the phase map is calculated by applying fast Fourier transform (FFT) to the time axis that is the frame direction for each pixel. However, maintaining real-time CD SWI is clinically important; real-time visualization of shear wave propagation on shear wave phase and velocity maps is required to allow physicians to properly position the vibrator and receiver device for image acquisition.

In this study, we investigated shear wave phase reconstruction using a limited number of CFIs. The proposed CD SWI phase map reconstruction method has two steps: the first step is to reconstruct an initial phase map from a limited number of CFIs using the advanced iterative algorithm (AIA)¹² used in optical interferometry; the second step is to reduce phase variations using the phase retrieval algorithm, based on an iterative procedure that cycles between Fourier and inverse Fourier space, while imposing a priori information on individual spaces.¹³ A directional filter (DF),^14-16 which extracts the shear wave components that propagate in the required direction, and total variation (TV) regularization¹⁷ are applied as constraints in individual spaces to suppress phase variations on the phase maps. We refer to this method as the directional filter and total variation (DFTV) method.

Bioucas-Dias et al. have proposed several algorithms related to de-noising or phase unwrapping of a phase map.^18,19 The PEARLS algorithm¹⁸ is based on local polynomial approximations, applying an adaptive filter while varying the window size in every pixel on a phase map and shows the excellent de-noising effect, even when the signal-to-noise ratio (SNR) is −0.5 dB. The PUMA algorithm¹⁹ is based on max-flow/min-cut calculations using the graph cut method to solve the integer optimization problem of phase unwrapping and can be reconstructed as an accurate unwrapped phase map for a noisy wrapped phase map.

Method

Acquisition of CFIs Including the Shear Wave Wavefront

A schematic diagram of the CD SWI system is shown in Figure 1. The system comprises a commercial ultrasound imaging system (EUB-8500 with a 6.5 MHz linear probe, Hitachi, Tokyo, Japan), a vibrator, and a personal computer (PC) for phase map reconstruction. The ultrasound probe and vibrator are positioned on the surface of the sample to be observed. When the vibrator, with a frequency that satisfies the CD SWI shear wave frequency condition, is placed on the sample surface, the generated shear wave propagates through the sample. The wavefront can be imaged using the CFI mode of an ultrasound imaging system. Captured CFI frames are immediately transferred to the PC through the ultrasound scanner video output for reconstruction. Using our proposed method, a single phase map is reconstructed from four CFI frames saved on the PC and is displayed on the screen. The reconstructed phase map is updated with the addition of each frame. For example, if the ultrasound imaging system CFI frames-per-second (FPS) is 10, the displayed phase map is refreshed every 100 ms. CFIs with a phase-shift of 0, 1/2π, π, and 3/2π radians are selected for phase map reconstruction because the phase map accuracy using the AIA method tends to be higher when the phase-shift interval between CFIs is equal within a wavelength.¹² The phase-shift interval is adjusted by shifting the shear wave frequency from the CD SWI frequency condition.

Figure 1.

Schematic diagram of the color Doppler shear wave imaging system.

Phase Map Reconstruction

When shear wave excitation leads to shear wave propagation inside the sample, the tissue moves sinusoidally by shear wave propagation. Subsequently, pseudo flow velocity patterns appear on the CFI because the high-frequency Doppler signal that is produced by shear wave propagation changes into a low-frequency signal via the aliasing process of sampling theorem. The pseudo flow velocity pattern represents the shear wave wavefront, which consists of the zero and maximum flow velocities on the CFI.¹¹ The brightness of the CFI showing the shear wave wavefront can be assumed as

v (x, y, i) = v_{\max} sgn (\cos [ϕ (x, y) + δ (i Δ t)]) where sgn (p) = {\begin{matrix} 1 & f o r & p \geq 0 \\ 0 & f o r & p < 0 \end{matrix},

where x and y are the pixel index in the lateral and axial directions, respectively; i (1 ≤ i ≤ N) and Δt are the CFI frame index and the time interval between neighboring frames, respectively; $ϕ$ and $δ$ are the phase in space and the amount of phase-shift in the frame direction, respectively; $v_{\max}$ is the maximum flow velocity on the CFI; and $v (x, y, i)$ includes the extensive artifacts from standing waves, reflected waves at tissue boundaries, and electric noise of the ultrasonic signal. Because Equation (1) is a periodic function of the pulse wave, the real and imaginary parts of φ(x, y) can be obtained by analyzing a fundamental wave component of $v (x, y, i)$ using FFT or AIA. Figure 2 shows a flowchart of the phase map reconstruction method:

Figure 2.

Flow chart of the shear wave phase map reconstruction method. AIA = advanced iterative algorithm; 2DFFT = two-dimensional fast Fourier transform.

Step 0. Prepare an initial signal map $f^{(0)} (x, y) = \cos ϕ (x, y) + j \sin ϕ (x, y)$ in which the real and imaginary parts are reconstructed from four CFIs using the AIA method.¹²

Step 1. Feed the initial signal map $f^{(0)} (x, y) = \cos ϕ (x, y) + j \sin ϕ (x, y)$ into the iterative procedure and set n = 1, where n shows the repetition number from Steps 1 to 6.

Step 2. Transform the updated signal map n times, $f^{(n)} (x, y)$ to the 2D Fourier map, $F^{(n)} (ξ, η)$ by applying the 2D FFT algorithm.

Step 3. Apply the DF to the 2D Fourier map. The DF, which has recently been used in other SWEI studies,^14
-16 separates the original forward-propagated shear wave from the other wave components by filtering the wave number vector map derived from the 2D Fourier transformation of the shear wave phase map. This reduces the phase artifacts caused by standing or reflected waves. However, a drawback of this filtering method is that it can be difficult to design a filter that only an original forward propagating shear wave can pass through; the propagation of shear waves through space, with refraction and diffraction at stiff tissue boundaries, is highly complex. Therefore, we designed a DF to extract shear waves that propagate away from the vibrator as a constraint condition in Fourier space. Thus, in the region S, that is, the right side of the 2D Fourier map, the previous pixel value is substituted for a new pixel, and in the region $\bar{S}$ , that is, the left side of the 2D Fourier map, zero is substituted for a new pixel.

{\tilde{F}}^{(n)} (ξ, η) = {\begin{matrix} F^{(n)} (ξ, η) & for & (ξ, η) \in S \\ 0 & for & (ξ, η) \in \bar{S} \end{matrix} .

Step 4. Transform the updated 2D Fourier space, ${\tilde{F}}^{(n)} (ξ, η)$ to the signal map, ${\tilde{f}}^{(n + 1)} (x, y)$ , using the inverse 2D FFT algorithm.

Step 5. Up to this step, the phase map contains phase variation because of the large region of $\bar{S}$ in Step 3; Step 5 suppresses the phase variations using a de-noising process based on TV regularization.¹⁷ This is one of the most popular and successful methodologies in image processing and is extensively used for processes such as image restoration and image deconvolution because of its excellent edge preservation. TV regularization is effective for suppressing phase variation because the phase-on-phase maps tend to change monotonically in the propagation direction. TV regularization is applied as follows:

f^{(n + 1)} (x, y) = \underset{f}{\arg \min} {‖ f - {\tilde{f}}^{(n + 1)} ‖}^{2} + λ TV (f),

where f is the desired true solution, λ controls the relative weights of the data fidelity and regularization terms, and

TV (f) = \iint_{X, Y} \sqrt{{(\frac{\partial f}{\partial x})}^{2} + {(\frac{\partial f}{\partial y})}^{2}} d x d y,

where X and Y are the size of the region of interest (ROI) in the lateral and longitudinal directions of the phase map, respectively. These sizes are similar to those of CFI images. The measurer can select the ROI of the CFI from the screen of the ultrasonic device. In this step, the objective functions of Equations (3) and (4) are minimized using the standard steepest descendent method and a derivative of the TV term. The step size and repetition number of the steepest descended method are important parameters in the DFTV method. In the simulation experiment, these values are set to n = n + 1.

Step 6. Iterate Steps 2 through 5 until

\frac{‖ f^{(n + 1)} (x, y) - f^{(n)} (x, y) ‖}{‖ f^{(n)} (x, y) ‖} \leq 10^{- 3} .

Finally, the phase map $ϕ (x, y)$ can be determined from

ϕ (x, y) = \tan^{- 1} [\frac{Im (f^{(n)} (x, y))}{Re (f^{(n)} (x, y))}],

where Re and Im denote the real and the imaginary parts of $f^{(n)} (x, y)$ respectively.

Experiments and Results

Simulation

To demonstrate the efficacy of the DFTV method, we perform a numerical simulation of CD SWI. A numerical phantom is a homogeneous medium with shear wave velocity of 3.0 m/s at the shear wave frequency of 275 Hz. The planar shear wave with shear wave amplitude of 65 µm and shear wave frequency of 275 Hz propagates from left to right into the phantom, as shown in Figure 3(a). Subsequently, the reflected planar shear wave with shear wave amplitude of 8 µm and shear wave frequency of 275 Hz propagates from right to left into the phantom. The ultrasonic probe for obtaining the CFI images is fixed on the top of the phantom. In the synthetic CFI image, the flow velocity on CFI is estimated from the quadrature detector output signals, which are acquired for M + 1 successive ultrasonic waves as²⁰

Figure 3.

(a) Diagram of the shear wave numerical simulation; (b) a CFI image used in the simulation experiment. CFI = color flow image; US = ultrasonic.

v = \frac{c}{4 π f_{0} Δ t} \arctan (\frac{E_{U}}{E_{L}}),

where

E_{U} = \sum_{m = 1}^{M} I_{m} Q_{m + 1} - I_{m + 1} Q_{m} and E_{L} = \sum_{m = 1}^{M} I_{m + 1} I_{m} + Q_{m + 1} Q_{m} .

f₀ is the center frequency of the ultrasonic wave, c is the sound velocity, and Δt is the time duration between successive ultrasonic pulses that are radiated in the same direction. In the simulation, f₀, c, Δt, M are set to 6.5 MHz, 1500 m/s, 2.74 ms, and 8, respectively. The I and Q signals have Gaussian noise with an SNR of 0 dB. The attenuation of the shear wave in the phantom is neglected. Figure 3(b) is a CFI image acquired by the simulation experiment. The interval of the shear wave wavefront is different spatially due to the standing wave generated by the propagation of the reflected wave. The matrix size of the CFI image is 128 × 128; the pixel size is 0.13 × 0.13 mm².

Figure 4(a) shows a theoretical phase map, and Figure 4(b) to (d) show phase maps reconstructed from four synthetic CFI images using FFT, AIA, and DFTV. The phase of each phase map is wrapped every π (radian) phase change. Figure 4(b) and (c) show strong phase variation, and the phase intervals are not equal to the phase modulation due to the standing wave on the CFI image. In contrast, Figure 4(d) shows strongly reduced phase variation due to the de-noising step, and the phase intervals are effectively corrected by the DF step method.

Figure 4.

Shear wave phase maps reconstructed from synthetic four CFI images. (a) Theoretical image, (b) FFT method, (c) AIA method, (d) DFTV method. CFI = color flow image; FFT = fast Fourier transform; AIA = advanced iterative algorithm; DFTV = directional filter and total variation.

We quantitatively evaluate the quality of the reconstructed phase maps. Figure 5(a) and (b) show the dependence of the root mean square error (RMSE) and the structural similarity (SSIM)²¹ on the number of CFIs. SSIM reaches a maximum value of 1.0 when the theoretical image and reconstructed image are entirely consistent. The RMSEs and SSIMs of FFT and AIA are worse than those of DFTV for every number of CFIs and begin to decrease at approximately 16 CFIs. In contrast, the RMSE and SSIM of DFTV are stable for four CFIs. Therefore, the quality of the DFTV is superior to those of FFT and AIA.

Figure 5.

Relationship between the image quality of a phase map and the number of CFI frames: (a) root mean square error, (b) structural similarity. CFI = color flow image; FFT = fast Fourier transform; AIA = advanced iterative algorithm; DFTV = directional filter and total variation.

The image quality and the speed of convergence in DFTV depend on the step size and repetition number of the steepest descendent method. It is difficult to determine the optimal parameters for every data set because the SNR of the CFI varies by imaging target, but parameters close to the optimal values can be determined by simulation. The repetition number of the steepest descendent method is fixed to 10 because early convergence of DFTV is required for real-time measurement. The remaining step size is determined based on the image quality of the reconstructed phase map. Figure 6(a) to (d) show the phase maps for step sizes of 10⁻¹, 10⁻², 10⁻³, and 10⁻⁴. Figure 6(a) and (b) are divergent due to an excessively large step size, and Figure 6(d) does not sufficiently remove phase variations due to a shortage of repetition. Figure 6(c) appears to be reconstructed properly. The RMSEs of Figure 6(a) to (d) are 1.349, 0.508, 0.123, and 0.284 radians, respectively. Therefore, we use a step size of 10⁻³ in all experiments.

Figure 6.

Reconstructed images for varying step size of the steepest descendent method in the proposed method: (a) 1.0 ×10⁻¹, (b) 1.0 ×10⁻², (c) 1.0 ×10⁻³, (d) 1.0 ×10⁻⁴.

Figure 7 shows the relationship between the quality of a phase map using DFTV and the SNR for the CFIs. The RMSE and SSIM abruptly begin to decrease at approximately 0 dB, which means that the noise exceeds the ultrasonic signal intensity. Based on this result, when the attenuation of the signal intensity is large, such as for adipose tissue, the quality of the phase map may degrade.

Figure 7.

Relationship between the image quality and SNR of the CFI images. SNR = signal-to-noise ratio; CFI = color flow image; SSIM = structural similarity; RMSE = root mean square error.

Experiments

Next, we demonstrate the efficacy of the DFTV method using actual experimental data of a commercial breast phantom and a human in vivo breast. The phantom, shown on the left side of the photograph in Figure 8, was designed for breast elastography training (OST Co., Ltd., BB-4, Ibaraki, Japan). The phantom contains models of adipose tissue, breast tissue with a cyst, and muscle tissue as shown by the right B-mode image in Figure 8. The purpose of the phantom is to reconstruct a phase map in tissues, except a cyst region. The human in vivo breast measurement was performed on a healthy 54-year-old woman. All experiments were approved by the institutional review board of Nissin Hospital.

Figure 8.

Photograph and B-mode image of the breast phantom.

A linear ultrasonic probe with a central frequency of 6.5 MHz was used in both reconstruction experiments. The shear wave was produced using a custom-made vibrator with a vibration amplitude of approximately 600 µm under an applied voltage of approximately 1 V. The resultant shear wave frequency was 276.5 Hz, which satisfies k = 1 of the shear wave frequency condition.¹⁰ The vibrator head and the ultrasonic probe were aligned to enable perpendicular propagation of the shear wave wavefront.

Results of the Breast Phantom Study

Figure 9 shows the four consecutive CFIs recorded using the breast phantom; these images were obtained from the area within the yellow square on the right side of Figure 8. The shear wave wavefront, propagating from left to right, is shown in red. The shear wave barely propagated in the phantom cyst, which is principally composed of liquid. The CFIs were processed using FFT, AIA, and DFTV. The phase map produced using AIA is similar to the output from Step 0 of the DFTV method. Figure 10(a) shows a phase map reconstructed using DFTV and 32 CFIs captured over approximately 3 s; this phase map is used as the standard for comparing the other phase maps reconstructed using four CFIs. Figure 10(b) to (d) show phase maps reconstructed from the four CFIs shown in Figure 4 using FFT, AIA, and DFTV. The shear wave did not propagate in the region of the cyst, which is colored black for easy recognition.

Figure 9.

The four CFIs acquired from the ultrasound device. CFI = color flow image.

Figure 10.

Shear wave phase maps reconstructed from the breast phantom CFIs: (a) DFTV method using 32 CFIs, (b) FFT method using four CFIs, (c) AIA method using four CFIs, (d) DFTV method using four CFIs. The dotted circle in the center indicates the cyst. CFI = color flow image; FFT = fast Fourier transform; AIA = advanced iterative algorithm; DFTV = directional filter and total variation.

Although phase changes should monotonically increase further from the vibrator, the phase maps shown in Figure 10(b), reconstructed using FFT, show a stepwise pattern. This pattern is not seen in Figure 10(c) or (d); we address this result in the “Discussion” section. In contrast to the phase variations observed in Figure 10(c) and (d), Figure 10(d) shows strongly reduced phase variations, which were achieved through the DFTV de-noising step. Figure 10(c), which shows the phase map reconstructed using AIA, shows rough structures in the vicinity of the phase wrap.

To evaluate the similarity between Figure 10(a) and Figure 10(b) to (d), we compared the line profiles obtained with each method, which are shown in Figure 11. The left and right columns show the profiles for the solid and dashed lines shown in Figure 10(a), respectively. The line profile produced using DFTV in Figure 11(c) is consistent with the profile of Figure 10(a); FFT and AIA, shown in Figure 11(a) and (b), are less consistent.

Figure 11.

Comparison of the line profiles from the different methods. The left and right graphs compare the line profiles from the solid and dotted lines in Figure 5(a), respectively. The individual graphs compare the line profile from Figure 5(a), obtained using the DFTV method and 32 CFIs (blue line), with the line profile obtained using four CFIs and (a) FFT, (b) AIA, or (c) DFTV (orange line). CFI = color flow image; FFT = fast Fourier transform; AIA = advanced iterative algorithm; DFTV = directional filter and total variation.

The RMSE and SSIM were calculated for all regions, excluding the cyst. The RMSEs of FFT, AIA, and DFTV were 1.27, 0.51, and 0.27 radians, respectively; the SSIMs of FFT, AIA, and DFTV were 0.02, 0.49, and 0.73, respectively.

Results of the In Vivo Human Breast Study

The shear wave observation region, enclosed in the red square of the B-mode image in Figure 12, includes adipose tissue (surrounded by an orange dotted line) and breast tissue. To reconstruct the phase map, four continuous CFIs were captured using the CD SWI system. Figure 13(a) shows the phase map reconstructed using DFTV and 32 CFIs captured over approximately 3 s. This map is used as the standard map for evaluating the other phase maps reconstructed from four CFIs. Figure 13(b) to (d) show the phase maps reconstructed from four CFIs using FFT, AIA, and DFTV. Figure 13(d), reconstructed by DFTV, shows the smallest phase variation. The region surrounded by the orange dashed line in Figure 13(a) to (d) and Figure 12 shows adipose tissue. In Figure 13(a) and (d), the adipose tissue region shows more complex shear wave propagations than the breast tissue. This might be a result of the non-uniform nature of adipose tissue, which includes complex fibrous components, or the low ultrasonic signal in the adipose tissue. In Figure 13(b) and (c), the adipose tissue region is unclear because of extensive noise. We calculated the phase velocity in the adipose tissue and breast tissue regions labeled A and B in Figure 13(b) to (d). These phase velocity values should be positive because the shear wave propagates from left to right on the map. The phase velocity of the 10 pixels adjacent to points A and B, corresponding to approximately 1.3 mm, are −4.06 m/s and −7.32 m/s in Figure 13(b), −3.68 m/s and −6.54 m/s in Figure 13(c), and 3.14 m/s and 2.76 m/s in Figure 13(d). The phase velocity values from Figure 13(b) and (c) are negative, indicating that large local phase variation prevented accurate phase velocity calculations. Only the phase velocity of Figure 13(d), reconstructed using DFTV, is positive. The shear wave velocities calculated for Figure 13(d) are within the range reported for adipose and benign breast tissue using ARFI: 3.05 ± 1.26 m/s and 5.39 ± 2.95 m/s, respectively.²² The shear wave velocity by Wojcinski et al is obtained from 88 cases and shows that the stiffness of the breast tissue has large variation due to individual differences. Because the shear wave velocity in Figure 13(d) is close to the lower value by Wojcinski et al, it is assumed that the human breast in the experiment is relatively soft. However, we should perform cross calibration between CD SWI and ARFI in the future. We also quantitatively evaluated the similarities between Figure 13(a) and Figure 13(b) to (d) by calculating the RMSEs and SSIMs of the entire images. The RMSEs of FFT, AIA, and DFTV were 1.69, 1.70, and 0.76 radians, respectively; the SSIMs of FFT, AIA, and DFTV were 0.12, 0.12, and 0.67. The RMSE values were generally quite high, indicating a low level of similarity; we assumed this to be attributed to respiration-related movement during the 3-s acquisition time required for Figure 13(a). Regardless, the RMSE and SSIM analysis indicates that DFTV is superior to FFT and AIA.

Figure 12.

B-mode image of the human breast. CFIs were acquired from the region of interest surrounded by the red square. CFIs = color flow images.

Figure 13.

Shear wave phase maps reconstructed from the human breast CFIs: (a) DFTV method using 32 CFIs, (b) FFT method using four CFIs, (c) AIA method using four CFIs, (d) DFTV method using four CFIs. CFIs = color flow images; FFT = fast Fourier transform; AIA = advanced iterative algorithm; DFTV = directional filter and total variation.

Discussion

We have shown that DFTV can retain image quality and reduce the number of CFIs required in the simulation experiment to four. DFTV was also effective in both phantom and human breast experiments. In this section, we discuss how the DFTV method using four CFIs speeds up reconstruction for real-time measurement, the benefits of using the AIA method as the initial phase map reconstruction method of DFTV, and the improvements gained by applying the DF constraint to each iteration of DFTV.

To evaluate the real-time imaging potential, we measured the time required to generate a single phase map, from the acquisition of four CFIs to the reconstruction of the phase map. The time taken to acquire four CFIs is approximately 200 to 266 ms using a routine ultrasound imaging system with a maximum 15 to 20 CFI FPS. The time required to reconstruct a phase map depends on the number of iterations and the computing power of the PC. Although the number of iterations may fluctuate with image size or the amount of CFI artifacts, we estimate the number of iterations from the convergence of the DFTV method in the breast phantom experiment. Figure 14 shows the plots of the left side of Equation (5) against the number of iterations. DFTV exhibits good convergence properties; using the criteria of Equation (5), the DFTV method converges to values less than 10⁻³ with seven or more iterations. The processing time required for seven iterations, including the time required for the initial phase map, is approximately 50 ms using a standard PC with a graphics processing unit. In total, the time required for CFI acquisition is approximately 300 ms; this is approximately one order of magnitude less than in previous studies using 32 CFIs.^10,11 This reduction could enable physicians to observe shear wave propagation on phase maps in near real-time.

Figure 14.

Convergence of the DFTV method. DFTV = directional filter and total variation.

DFTV employs AIA for initial phase map reconstruction. AIA reconstructs accurate phase maps from at least three frames with different wavefronts by repeating least-squares calculations between the time and space domains. If the time-domain phase-shift interval is different, the phase map can be estimated using the space-domain wavefront information. In the CFIs of Figure 4, the phase-shift interval of the four CFIs is slightly shifted from 0, ½π, π, and 3/2π radians because of the different FPS in the video output and the CFIs. Therefore, the time-domain fundamental frequency was different from the shear wave frequency. This is illustrated by the stepwise artifacts in Figure 10(b), which was generated using FFT. Conversely, Figure 10(c) and (d), obtained using AIA, do not show stepwise artifacts because both time and space information was used. Therefore, AIA is more effective than FFT for initial phase map reconstruction.

The noise removal method based on TV minimization can extensively remove noise from the image. However, the structure of the image is over-flattened. When applying it to a phase map with monotonously increasing values, a false image with a stepwise pattern is generated. Figure 15(a) shows a phase map with TV minimization applied to the phase map of Figure 4(c) in the simulation section. For comparison, the step size and the repetition number of the steepest descended method are similar to the conditions used for Figure 4(d). The phase variation in Figure 15(a) is decreasing, while the phase distribution does not monotonously increase but shows stepwise shapes. In contrast, in the Fourier space of this phase map (i.e., the wave number vector map), the stepwise shape component is covered in a wide bandwidth because a stepwise shape is represented as an infinite Fourier series. Since the DF removes the component of reverse propagation (i.e., one side of the wave number vector map), stepwise shapes on the phase map can be corrected to some extent, even if the DF is applied only once to Figure 15(a), as shown in Figure 15(b). However, Figure 15(b) has more noise than Figure 4(d). Figure 15(c) shows a phase map–applied TV minimization after applying the DF to Figure 4(c). The irregular phase interval due to the standing wave is corrected, but stepwise shapes are observed. The RMSEs of Figure 15(a) to (c) are 0.387, 0.154, and 0.171 radians, respectively, and the SSIMs are 0.547, 0.846, and 0.829, respectively. Because the RMSE and SSIM of Figure 4(d) are 0.123 radians and 0.840, respectively, the RMSEs of Figure 15(b) and (c) are slightly larger than that of Figure 4(d). Therefore, the application of the DF to each iteration of DFTV is effective to correct the adverse influence of TV minimization.

Figure 15.

(a) Phase map applying the total variation minimization. (b) Phase map applying a directional filter to the phase map of (a). (c) Phase map applying the total variation minimization after applying a directional filter to the initial phase map.

Some limitations of the current study should be highlighted. First, our study was conducted with a small number of simple structure imaging targets; the spatial resolution and applicability of the DFTV method need to be verified using phantom or biomedical models with more complex structures. In addition, the accuracy of the phase map reconstructed using the DFTV four-CFIs method was demonstrated by comparing it with a 32 CFI phase map. This assumes that the 32 CFI phase map is correct. To verify the reference phase map, simulation experiments using numerical phantoms and simulation data, including standing waves or waves reflected from a tissue boundary, must be conducted.

Conclusion

We have proposed a method for phase map reconstruction using a limited number of CFIs acquired with a CD SWI system. We have demonstrated the efficacy of this method using synthetic and experimental data of a breast phantom and human breast tissue. Our results show that the proposed DFTV method can maintain image quality and reduce the number of CFIs required to four; previous methods have required at least 32 CFIs to achieve equivalent image quality. The DFTV method is suitable for real-time SWEI using a continuous shear wave produced by a mechanical vibrator.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was supported by a Grant-In-Aid for Scientific Research (B) #16H04373 from the Ministry of Education, Culture, Sports, Science and Technology in Japan.

References

Castera

Foms

Alberti

Non-invasive evaluation of liver fibrosis using transient elastography. J Hepatol. 2008;48:835-47.

Itoh

Ueno

Tohno

Kamma

Takahashi

Shiina

. Breast disease: clinical application of US elastography for diagnosis. Radiology. 2006;239:341-50.

Tohno

Ueno

Current improvements in breast ultrasound, with a special focus on elastography. Breast Cancer. 2008;15:200-4.

Hong

Liu

Zhang

Chen

Luo

Real-time ultrasound elastography in the differential diagnosis of benign and malignant thyroid nodules. J Ultrasound Med. 2009;28:861-7.

O’Donnell

Skovoroda

Shapo

Emelianov

SY.

Internal displacement and strain imaging using ultrasonic speckle tracking. IEEE Trans Ultrason Ferroelectr Freq Control. 1994;41:314-25.

Nightingale

Soo

Nightingale

Trahey

Acoustic radiation force impulse imaging: in vivo demonstration of clinical feasibility. Ultrasound Med Biol. 2002;28:227-35.

Bercoff

Tanter

Fink

Supersonic shear imaging: a new technique for soft tissue elasticity mapping. IEEE Trans Ultrason Ferroelectr Freq Control. 2004;51:396-409.

Yamakawa

Nitta

Shiina

Matsumura

Tamano

Mitake

. High-speed freehand tissue elasticity imaging for breast diagnosis. Jpn J Appl Phys. 2003;42:3265-70.

Nagaoka

Iwasaki

Arakawa

Kobayashi

Yoshizawa

Umemura

. Basic study of intrinsic elastography: relationship between tissue stiffness and propagation velocity of deformation induced by pulsatile flow. Jpn J Appl Phys. 2015;54:07HF08.

10.

Yamakoshi

Kasahara

Iijima

Yuminaka

Shear wave wavefront mapping using ultrasound color flow imaging. Ultrason Imaging. 2015;37:323-40.

11.

Yamakoshi

Yamamoto

Kasahara

Iijima

Yuminaka

Shear wave mapping of skeletal muscle using shear wave wavefront reconstruction based on ultrasound color flow imaging. Jpn J Appl Phys. 2015;54:07HC16.

12.

Wang

Han

Advanced iterative algorithm for phase extraction of randomly phase-shifted interferograms. Opt Lett. 2004;29:1671-3.

13.

Fienup

JR.

Reconstruction of an object from the modulus of its Fourier transform. Opt Lett. 1978;3:27-29.

14.

Deffieux

Gennisson

Bercoff

Tanter

On the effects of reflected waves in transient shear wave elastography. IEEE Trans Ultrason Ferroelectr Freq Control. 2011;58:2032-5.

15.

Ono

A method to reduce the influence of reflected waves on shear velocity measurements using B-mode scanning time delay. Jpn J Appl Phys. 2015;54:07HF01.

16.

Sarvazyan

Urban

Greenleaf

JF.

Acoustic waves in medical imaging and diagnostics. Ultrasound Med Biol. 2013;39:1133-46.

17.

Rudin

Osher

Fatemi

Nonlinear total variation based noise removal algorithms. Physica D. 1992;60:259-68.

18.

Bioucas-Dias

Katkovnik

Astola

Egiazarian

Absolute phase estimation: adaptive local denoising and global unwrapping. Appl Opt. 2008;47:5358-69.

19.

Bioucas-Dias

Valadão

Phase unwrapping via graph cuts. IEEE Trans Image Process. 2007;16:698-709.

20.

Kasai

Namekawa

Koyano

Omoto

Real-time two-dimensional blood flow imaging using an autocorrelation technique. IEEE Trans Ultrason Ferroelectr Freq Control. 1985;32:458-64.

21.

Wang

Bovik

Sheikh

Simoncelli

EP.

Image quality assessment: from error visibility to structural similarity. IEEE Trans Image Process. 2004;13:600-12.

22.

Wojcinski

Brandhorst

Sadigh

Hillemanns

Degenhardt

Acoustic radiation force impulse imaging with Virtual Touch™ tissue quantification: mean shear wave velocity of malignant and benign breast masses. Int J Womens Health. 2013;5:619-27.