Two-dimensional complex shear modulus imaging of soft tissues by integration of Algebraic Helmoltz Inversion and LMS filter into dealing with noisy data: a simulation study

1.   Introduction

Elasticity and viscosity are two tissue characteristic parameters; they can be used to detect tumors in the human body ^[1]. In the last few decades, a number of researchers have adopted a detection technique based on the tissue elasticity ^{[2,3,4,5,6,7]}, namely the ultrasonic Shear Wave Elasticity Imaging (SWEI). SWEI offers the possibility to image the elasticity of tissues and organs that the traditional ultrasound cannot reach. Meanwhile, for a better understanding of the tissues, it would be beneficial to simultaneously look at these two parameters due to the fact that certain tissues have similar elasticity and different viscosity, for example, in Liver ^[1]. Many authors have developed different estimation techniques for both the elasticity and viscosity ^{[7,8,9,10,11,12,13,14,15,16]}.

Ultrasound Imaging and other methods such as Mammography, Clinical Examination and MR Imaging can be used to diagnose the breast cancer ^[17]; among these, shear wave has recently been employed to estimate CSM. In practice, the vibration force used to generate shear wave is in the range $50-250$ Hz. In ^[15], the authors conducted an experiment in which they estimated the elasticity and viscosity of a phantom sample. The true values of the elasticity and viscosity of a phantom sample were known. The results showed that the true and the estimated values of the elasticity and viscosity are quite closed at 150 Hz. In this paper, we propose using the same frequency of $150$ Hz in order to compare our results to the results of this study. The relationship between the speed and absorption of shear wave is modeled by Greanleaf in ^[18]. In order to minimize the reflections, several authors proposed using pulse excitation ^[19]. However, the influence of noise was more adverse than using harmonic ones ^[20]. Another advantage of using harmonics is that a much larger amplitudes of shear wave can be obtained. For these reasons, we used harmonic excitation in our work.

In 2004, the shear wave speed dispersion was established by Chen et al. using the relationship between the propagation speed and the vibrating frequency; this speed dispersion was finally used to estimate the CSM ^[8]. In 2007, Zheng et al. used the same technique proposed by Chen but added a linear Kalman filter to minimize noise influence ^[20]. Recently, this method was used by several authors ^[8,10,15,20] where multiple datasets of distinct vibration frequencies must first be obtained. In this paper, only a single-frequency excitation was needed, which tremendously improved the acquisition time. Tran et al. ^[21] used the maximum likelihood ensemble filter for $1D$ to detect tumors (if any) in the tissues. However, the propagation model using wave equation in ^[10,14,21] is too simple to represent the actual propagation in heterogeneous tissues. In 2015, Qiuang et al. ^[22] proposed a technique using Finite Element Method (FEM) to better produce numerical solutions for the transversal shear wave propagation in any medium. However, FDTD is a more efficient method with lower memory usage than FEM ^[23]. In ^[12], Orescanin et al. combined the FDTD method with the AHI algorithm to estimate CSM. In this study, authors used a low-pass second-order Butterworth filter to reduce the noise impact.

In our previous study ^[24], we integrated Least Mean Square (LMS) with AHI algorithms to estimate 1D CSM. In fact, the 1D CSM estimation did not serve much purpose in the diagnosis. But the findings of this study confirmed that the noise reduction of acquired data before estimating CSM will produce better results. Thus, in this study, we proposed a method for 2D-CSM imaging of soft tissues based on the integration between AHI and Normalize Least Mean Square (NLMS) algorithms. The simulation scenarios were carried out in the MATLAB environment. An actuator was used to vibrate a needle at the frequency of $150$ Hz. This needle acted on the surface of tissue, leading to the generation and propagation of shear wave in the tissues. A Doppler device was then used to measure the particle velocity of shear wave at $120 \times 120$ spatial locations. At each spatial location, the noise coming from the measured particle velocity was reduced noise using NLMS. After that, CSM was directly estimated by using the AHI algorithm. Finally, the 2D-CSM image was reconstructed. Our proposed method NLMS/AHI offers a very low complexity compared to the previous methods. A scenario with a tumor and noise environment was studied to evaluate the quality of the estimated CSM.

2.   Method

2.1.   Simulation model

The measurement system consisted of a vibrating needle, an actuator, and a Doppler ultrasound device (see Figure 1). The actuator generated through a needle a mechanical vibration. Finally, we assumed that the Doppler ultrasound device was used to measure the shear wave at each spatial location in the region of interest ^[10]. In the later numerical simulation, the particle velocity which propagated in tissues was modeled by using FDTD model. This particle velocity, after adding noises, simulated the measured particle velocity (from the Doppler device). Thus the measured particle velocity was the noisy particle velocity.

Figure 1.  Generation and measurement of shear wave.

DownLoad: Full-Size Img PowerPoint

In this paper, we estimated the elasticity $\mu$ and the viscosity $\eta$ of the tissues, which were based on the Kelvin-Voigt model as follows

where $\omega$ is the angle frequency of the vibration signal, $\mu$ is the elasticity and $\eta$ is the viscosity at the spatial position (x, y).

A system of six first-order hyperbolic equations was used to describe wave propagation in 2D space as shown in the following equations (2.2), (2.3), and (2.4) ^[12].

where $v_z$ is the particle velocity vector on the plane (x, y), $\sigma _{zx}$ and $\sigma _{zy}$ are two elements of the stress tensor $\sigma$, $\partial _t$ denotes a partial derivative operator $\partial / \partial t$ which is applied to the parameters to the right of the symbol, $\partial _x$ denotes a partial derivative operator $\partial / \partial x$ applied to the parameters to the right of the symbol, $\partial _y$ denotes a partial derivative operator $\partial / \partial y$ applied to the parameters to the right of the symbol, $\rho$ is the density of the tissues.

We discretized $v_z$, $\sigma _{zx}$ and $\sigma _{zy}$ in Eqs (2.2), (2.3) and (2.4) to obtain the following equations

where $\Delta x$ and $\Delta y$ are the distances between continuous spatial locations following the X-axis, and Y-axis, respectively; $\Delta t$ is the sampling cycle, the index $i$ denotes the spatial step following the X-axis, the index $j$ denotes the spatial step following the Y-axis, and the index $n$ denotes the temporal step, These parameters are illustrated in Figure 2.

Figure 2.  The particle velocity and stress tensor nodes on plane–(X, Y).

DownLoad: Full-Size Img PowerPoint

By applying the FDTD method in 2D medium, Equations (2.2), (2.3) and (2.4) are transformed into Equations (2.8), (2.9) and (2.10), respectively.

The accuracy of numerical Eqs (2.8), (2.9) and (2.10) is the first order in time and space. The boundary conditions were applied by using an absorption boundary layer. It will reduce the reflections. The layer is defined by a loss profile chosen within a perfectly matched absorptive layer and the attenuation factor $\tau$ which is given by ^[25]

where two constants are set as $m = 2.15$ and $\tau _{m{\rm{ax}}} = 0.11$; the layer thickness in grid points $N_{pml} = 10$; and the term $i - i_{pml}$ indicates position within the perfectly matched absorptive layer.

2.2.   Signal improvement using our proposed approach

FDTD and AHI can be used to directly estimate CSM at each point in the space (presented in detail in section 2.3). In such an approach, one important point is to minimize the propagation of error, which mean that the influence of measurement noise need to be reduced as much as before the direct estimation. In our work, we propose to use adaptive filtering ^[26] which was proved to give a mean error of $2.06 \%$.

Figure 3 describes the LMS adaptive filter we used to eliminate noise from the noisy signal. The adaptive filtering algorithm is described in Algorithm 1. The noisy particle velocity generation can be described as follows: the Eqs (2.8), (2.9) and (2.10) were used to create the ideal particle velocity $d(n)$ in each spatial location. This signal was added by a Gaussian white noise $v(n)$ to form a noisy particle velocity ${v_z}(n)$. The level of noise which affected to the particle velocities at the different observed locations was identical. Thus, the signal to noise ratio (SNR) reduced from the locations near to the needle to the ones farthest away from the needle. In our simulation scenario, the average SNR of all the spatial locations was $10.1$ dB. The reference noise generation can be described as follows: when the vibration needle stayed still, the measured signal was just the reference noise $x(n)$ as indicated in Figure 3. Note that $x(n)$ was the Gaussian white noise that we measured when the vibration needle stayed still, and $v(n)$ was the Gaussian white noise that we measured when the vibration needle was active. In Figure 3, $Var$ variable was the estimated variance of noise in ${v_z}(n)$. Both $x(n)$ and ${v_z}(n)$ were applied to a filter in ^[27] (chapter 9, pages 505 to 515).

Figure 3.  Signal improvement using NLMS filter.

DownLoad: Full-Size Img PowerPoint

Algorithm 1: NLMS algorithm for signal improvement

Step 1: Set up the parameters: step size $\beta$, filter order $L$, and noise variance. Step 2: Start the simulation by assigning $w(n) = 0$. Step 3: For $n = 0, 1, 2, \ldots$. $3.1.$ Calculate the output $y(n)$ $y(n) = w(n)*x(n)$ $3.2.$ Calculate the error $e(n)$ $e(n)$ = ${v_z}(n)$ - $y(n)$ $3.3.$ Update the filter's coefficient $w(n+1)$ = $w(n)$ + $\beta$ $e(n)x(n)$ $3.4.$ Extract the improved signal ${\widehat v_z}(n)$ ${\widehat v_z}(n)$ = $e(n)$ $3.5.$ Change the value of $\beta$ based on the amplitude of ${v_z}(n)$

---

 | Show Table

DownLoad: CSV

In practice, the shear wave was quickly attenuated with the distance from the position of the vibration needle. Due to this fact, the Signal-to-Noise Ratio (SNR) varied with the space points where we would like to estimate the CSM. The optimal step size, giving a trade-off between the speed of convergence and residual error, depended on the SNR. Therefore, the commonly used modification of the LMS algorithm was a normalized version of the step size. This leads to the normalized LMS (NLMS) algorithm, in which the step size was scaled (divided) by an estimated power of the input data (i.e the variable $Var$ as shown in Figure 3). This observation required that we use a different step size $\beta$. The complexity of this process suggested that we should use a simple structure for adaptive filtering with fast convergence ^[28].

The NLMS algorithm which was used in this study is shown in Algorithm 1 below. We adjusted the step size at each spatial location using the formula $\beta = \frac{{\beta _1}}{{\left({\varepsilon + Var} \right)}}$, where $\beta_1$ and $\varepsilon$ are constant, and $Var$ is the variance of ${v_z}(n)$ which was estimated at each spatial location as shown in Figure 3.

2.3.   2D–CSM Imaging using NLMS/AHI algorithm

The filtered particle velocity was brought to the AHI algorithm ^[29] to directly estimate the CSM. We assumed that the viscoelastic properties of the medium were isotropic and there was negligible compression applied to the medium by the source. Thus, the particle velocity vector $v_z$ can be described by the Navier wave equation in a homogeneous solid ^[12]. By combining (2.2), (2.3) and (2.4):

where $G'\left({x, y, t} \right)$ is the CSM in the time domain and ${\nabla ^2}{v_z}$ is the Laplace operator of $v_z$ which is defined as ${\nabla ^2}{v_z} = {\partial ^2}{v_z}/\partial {x^2}$.

AHI algorithm was then used to solve (2.12), which then converted to the Helmholtz equation

where $G\left({x, y, \omega } \right)$ is the CSM in the frequency domain and is defined in (2.1), ${V_z}\left({x, y, \omega } \right)$ is the temporal Fourier transform of the particle velocity ${v_z}\left({x, y, t} \right)$, ${V_z}\left({x, y, \omega } \right) = {F_t}\left\{ {{v_z}\left({x, y, t} \right)} \right\}$, and $\omega _0$ is the angular frequency ${\omega _0} = 2\pi f_0$. Solving (2.13), we could directly estimate CSM

Following (2.1), $\mu \left(x, y \right)$ and $\eta \left(x, y \right)$ were real and imaginary parts of $G\left({x, y, \omega } \right)$. Thus, elasticity and viscosity of the tissues at each spatial point $(x, y)$ were calculated as

where ${V_z}({x, y, {\omega _0}})$ was computed by using Fourier transform at the specific angular frequency ${\omega _0}$; ${\nabla ^2}{V_z}({x, y, {\omega _0}})$ was computed by using the function Discrete Laplacian (The MathWorks) $del2({V_z}({x, y, {\omega _0}}))$ which provided a discrete approximation of Laplace's differential operator applied to ${V_z}({x, y, {\omega _0}})$.

Our proposed NLMS/AHI algorithm is summarized in Algorithm 2.

Algorithm 2: NLMS/AHI algorithm for 2D-CSM imaging

Step 1: Initiate the parameters for simulation. Step 2: Vibrate the needle at a frequency of ${f_0}$=$150$Hz. Step 3: Propagate the shear wave into the tissues. Step 4: Use the Doppler ultrasound device to measure $v_z$ at $120 \times 120$ spatial locations. Step 5: Extract the variance of noise and forward to the NLMS filter. Step 6: Employ NLMS filter (as shown in Algorithm 1) in order to improve the SNR. Step 7: Compute the temporal Fourier transform of the filtered signal. Step 8: Directly estimate CSM in each spatial location using (2.15). Step 9: Reconstruct the 2D-CSM image.

---

 | Show Table

DownLoad: CSV

3.   Simulation and results

3.1.   Simulation

The researchers set up a simulation scenario to test the proposed method. Details are as follows: The 2D medium had a size of $120 \times 120\; mm$, which contained a tumor at the coordinate $(40\; mm, 40\; mm)$, the radius of the tumor was $20\; mm$. The elasticity of the medium was $\mu_1 = 6000\; Pa$ and the viscosity of the medium was $\eta = 1.2\; Pa.s$, The elasticity and viscosity of the tumor were $\mu_1 = 9000\; Pa$ and $\eta = 1.8\; Pa.s$. The parameters of the medium and tumor are summarized in Table 1.

Table 1.  Parameters of medium and tumor.

	Coordinate (mm, mm)	Radius (mm)	$\mu_1$ (Pa)	$\eta$ (Pa.s)
Medium			6000	1.2
Tumor	(40, 40)	20	9000	1.8

---

 | Show Table

DownLoad: CSV

The mass density of the medium was $\rho = 1000\; Kg/m^3$. On the 2D plane of the medium (tissue), the particle velocity of shear wave was measured at locations which were equally spaced of $1\; mm$ following both the X–axis and Y–axis. Thus, the number of measured locations was $120 \times 120 = 14400$. Figure 4 illustrates the ideal–elasticity and ideal–viscosity images.

Figure 4.  Ideal 2D image for the elasticity and viscosity of tissue.

DownLoad: Full-Size Img PowerPoint

3.2.   Results

The ideal particle velocity at a location was modeled following the FDTD method (as shown in Figure 5 at the spatial location (56, 56)). Combining with the Gaussian noise, it presented the measured particle velocity using the Doppler system.

Figure 5.  The ideal particle velocity at the spatial location (56, 56) which is zoomed in the range of 0.48 to 0.52 second.

DownLoad: Full-Size Img PowerPoint

Figure 6 indicates two signals of the particle velocity of the shear wave at the coordinate (56, 56). where the dotted–blue signal shows the measured particle velocity (affected by noise), while the continuous–red signal shows the filtered particle velocity.

Figure 6.  The filtered particle velocity at the spatial location (56, 56) which is zoomed in the range of 0.48 to 0.52 second.

DownLoad: Full-Size Img PowerPoint

Figures 7(a) and 7(b) show estimated–elasticity and estimated–viscosity images of the tissues, respectively. In both images, tumor is seen clearly in in terms of location and shape.

Figure 7.  Estimated 2D image for the elasticity and viscosity of tissue.

DownLoad: Full-Size Img PowerPoint

Observing 1D–data is easy to compare between the estimated and ideal values. Figure 8 illustrates results of elasticity estimation on Line 10 and Line 40. In detail, Figure 8(b) indicates the result of elasticity estimation on Line 10 which does not across tumor, Figure 8(c) shows the result of elasticity estimation on Line 40 which crosses tumor.

Figure 8.  Results of elasticity estimation on Line 10 and Line 40.

DownLoad: Full-Size Img PowerPoint

Figure 9 illustrates the results of viscosity estimation on Line 10 and Line 40. In detail, Figure 9(b) indicates the result of viscosity estimation on Line 10, Figure 9(c) shows the result of viscosity estimation on Line 40.

Figure 9.  Results of viscosity estimation on Line 10 and Line 40.

DownLoad: Full-Size Img PowerPoint

Observe Figures 8(b), 8(c), 9(b) and 9(c), it is clearly seen that the results of estimation were less accurate at locations which were far from excited source (the right of figures). This can explain as follows as shear wave was attenuated when it propagated in the tissues and it was a damped oscillation. While CSM was estimated base on measurement of velocity of shear wave. At locations which were far from excited source, velocity value was very small, thus CSM estimation was less accurate.

3.3.   Performance

To evaluate the quality of 2D CSM estimation, the normalized error was used. This norm which was applied to elasticity and viscosity images was calculated using the following equation (3.1).

where $\epsilon_{\mu}$ and $\epsilon_{\eta}$ are the normalized errors of 2D elasticity and viscosity estimation, respectively, $M \times N$ is the size of the images, $\mu_{i, j}$ and $\hat{\mu}_{i, j}$ are the ideal and estimated elasticity at the pixel $(i, j)$, respectively, $\eta_{i, j}$ and $\hat{\eta}_{i, j}$ are the ideal and estimated viscosity at the pixel $(i, j)$, respectively.

The calculated results show that the normalized error $\epsilon_{\mu}$ for 2D elasticity estimation was $1.63 \%$, while the normalized error $\epsilon_{\eta}$ for 2D viscosity estimation was $2.49 \%$.

4.   Discussions

It is noteworthy that the estimated images contain many interesting features. In Figures 7(a) and 7(b), near the vibration needle, the reconstructed image is of very high quality, i.e. the reconstructed image looks like the original one (used in the simulation scenario). The area where the estimation gives good result is the circle of radius of $120$ mm centered around the needle position. At spatial locations outside of this area, the estimation is not accurate because of the high attenuation of the particle velocity. Combining different estimates of CSM with the varying positions of the vibration needle is a good subject for future research.

The propagation of noise in a scheme such as AHI and FDTD is a serious problem. Our research showed that by filtering noise with NLMS, the situation was nicely ameliorated and the estimation estimation error was reduced to $1.63 \%$ for elasticity and $2.49 \%$ for viscosity. In the existing research works, for instance ^[12], the noise impact was reduced by applying a low-pass filter. With the same simulation condition, in our study, the normalized error for 2D elasticity estimation was $4.15 \%$, and for 2D viscosity estimation was $6.24 \%$. The reason is that, the low pass filter used in ^[12] could only remove the high frequency noise.

While there are significant findings reported in our study, the study itself suffers from certain limitations. Firstly, due to the variation of SNR in the measurement at each point in the space, if the step size of the adaptive filtering is kept fixed, the efficiency of the filter would be reduced. In our work, this parameter was changed depending on the position of the needle. Secondly, certain sets of data obtained from the measurement need to be eliminated. In fact, to estimate CSM using AHI, it requires that the Fourier transform ${V_z}({x, y, {\omega _0}})$ be at the frequency $\omega _0$. The transient in the measurement process could badly affect the estimation of ${V_z}({x, y, {\omega _0}})$; consequently, the transient has to be eliminated before we estimate ${V_z}({x, y, {\omega _0}})$ (see step 7, Algorithm 2).

In addition, in the shear wave propagation model, we assumed that the reflections and the refractions from the surrounding tissues are trivial, thus we ignored them. Finally, this research is only concerned with numerical simulation. The estimation using experimental data is our next target.

5.   Conclusion

This paper presented a method of constructing 2D CSM images using the estimation method AHI in a noisy environment. In our method, we just measured the particle velocity at only one frequency ($150$ Hz) at different positions in the focused area. Noise in the measurement of particle velocity was eliminated by using a NLMS adaptive filter in order to improve the image reconstruction. The quality of image reconstruction was improved by (1) Changing the adaptive filter step size at each point of the observation space and ($2$) Eliminating the transient data before computing the Fourier transform of the particle velocity at the frequency $150$ Hz; this helped to improve the estimation of CSM using Eq (2.15). These estimates computed at each point in the space as illustrated in Figure 2. were finally used to reconstruct the CSM images. Our proposed method delivered a good performance by reducing the estimation error to $1.63 \%$ for elasticity and to $2.49 \%$ for viscosity. The next steps in our future research will be applying our proposed method to real data and extending it to 3D CSM imaging.

Acknowledgments

The authors would like to sincerely thank Professor Michael F. Insana from the University of Illinois at Urbana–Champaign for introducing us to shear wave imaging. This work was supported by the Asia Research Center (ARC), Vietnam National University, code CA.17.6A.

Conflict of interest

The authors declare there is no conflicts of interest.