MR Elastography with Optimization-Based Phase Unwrapping and Traveling
  Wave Expansion-based Neural Network (TWENN)
Magnetic Resonance Elastography (MRE) can characterize biomechanical
properties of soft tissue for disease diagnosis and treatment planning.
However, complicated wavefields acquired from MRE coupled with noise pose
challenges for accurate displacement extraction and modulus estimation. Here we
propose a pipeline for processing MRE images using optimization-based
displacement extraction and Traveling Wave Expansion-based Neural Network
(TWENN) modulus estimation. Phase unwrapping and displacement extraction were
achieved by optimization of an objective function with Dual Data Consistency
(Dual-DC). A complex-valued neural network using displacement covariance as
input has been constructed for the estimation of complex wavenumbers. A model
of traveling wave expansion is used to generate training datasets with
different levels of noise for the network. The complex shear modulus map is
obtained by a fusion of multifrequency and multidirectional data. Validation
using images of brain and liver simulation demonstrates the practical value of
the proposed pipeline, which can estimate the biomechanical properties with
minimum root-mean-square-errors compared with state-of-the-art methods.
Applications of the proposed method for processing MRE images of phantom,
brain, and liver show clear anatomical features and that the pipeline is robust
to noise and has a good generalization capability.
Shengyuan Ma, Runke Wang, Suhao Qiu, Ruokun Li, Qi Yue, Qingfang Sun, Liang Chen, Fuhua Yan, Guang-Zhong Yang, Yuan Feng
MR Elastography with Optimization-Based Phase Unwrapping and Traveling Wave Expansion-based Neural Network (TWENN)

###### Abstract

Magnetic Resonance Elastography (MRE) can characterize biomechanical properties of soft tissue for disease diagnosis and treatment planning. However, complicated wavefields acquired from MRE coupled with noise pose challenges for accurate displacement extraction and modulus estimation. Using optimization-based displacement extraction and Traveling Wave Expansion-based Neural Network (TWENN) modulus estimation, we propose a new pipeline for processing MRE images. An objective function with Dual Data Consistency (Dual-DC) has been used to ensure accurate phase unwrapping and displacement extraction. For the estimation of complex wavenumbers, a complex-valued neural network with displacement covariance as an input has been developed. A model of traveling wave expansion is used to generate training datasets for the network with varying levels of noise. The complex shear modulus map is obtained through fusion of multifrequency and multidirectional data. Validation using brain and liver simulation images demonstrates the practical value of the proposed pipeline, which can estimate the biomechanical properties with minimal root-mean-square errors when compared to state-of-the-art methods. Applications of the proposed method for processing MRE images of phantom, brain, and liver reveal clear anatomical features, robustness to noise, and good generalizability of the pipeline.

 Magnetic resonance elastography, Modulus estimation, Neural network, Traveling waves, Phase unwrapping

## 1 Introduction

Magnetic Resonance Elastography (MRE) can measure viscoelastic mechanical parameters of soft tissues noninvasively [1, 2]. In recent years, studies have investigated the diagnostic potential of MRE for liver cirrhosis [3], brain tumors [4], liver tumors [5], and Parkinson's diseases [6]. Post-processing of MRE images includes two major steps: 1) extraction of wavefield from the original wrapped phase; and 2) estimation of biomechanical properties based on the measured wavefield [7]. However, accurate displacement extraction from wrapped phase images with noise remains difficult [8]. In vivo wavefields of soft tissues with complicated anatomical structures usually contain deflections and reflections with acquisition noise. These bring challenges to the accurate estimation of the biomechanical properties of soft tissues [9].

As the first step of processing MRE images, the quality of the extracted wavefield determines the overall performance of the extracted biomechanical properties. In phase images, wave information is encoded with a motion encoding gradient, such that the magnitude of the motion can introduce a wrapping effect. Thus, phase unwrapping is necessary. Conventional approaches involve unwrapping a single image frame before performing Fast Fourier transform (FFT) in the time domain to extract the principal component [10]. In addition, spatiotemporal information can also be used for phase unwrapping [11]. Unwrapping algorithms such as Sorting by Reliability (SG) algorithm [12] and Dilate-Erode-Propagate (DE) [8] work in a search-like or discrete greedy update manner, which can fail in complicated wrapped scenarios [8]. Laplacian-based Estimation (LBE) simultaneously performs phase unwrapping and FFT in the frequency domain [13]. LBE is noise-resilient, but may introduce additional background offset noise. This can have a negative impact on the estimated biomechanical properties. The phase gradient method (PG) does not need direct phase unwrapping, but is noise-sensitive and requires specially designed modulus estimation algorithms [14]. Therefore, a noise-resistant, accurate method for wavefield extraction that can perform both phase unwrapping and extraction of the principal component is required.

MRE is a shear-wave-based elastography that requires the estimation of complex wavenumbers from acquired wave images. Numerous methods have been proposed to date [9] and
typical algorithms include Direct Inversion (DI) using the Helmholtz equation [15] and local frequency estimation (LFE) [16]. Due to the use of Laplace operators, DI is susceptible to noise and prone to edge artifacts. Iterative algorithms with prior information such as Multifrequency Elasticity Reconstruction using Structured Sparsity and ADMM (MERSA) [17] and MRE Inversion by Compressive Recovery (MICRo) [18] were proposed to use DI as an estimation kernel. But these methods had similar problems due to the differentiation operation. LFE is more noise-resistant, but its ability to recover anatomical features is limited. Enhanced Complex Local Frequency (EC-LFE) was developed to provide viscosity information, but structural resilience remained a problem [19].

To overcome the limitations of the DI-based method, Multifrequency Dual Elasto-Visco inversion (MDEV) was proposed to improve noise robustness using multi-frequency data averaging [20]. However, differentiation operation could still introduce structural edge artifacts, making it sensitive to noise. To improve the capability of distinguishing anatomical features, k-MDEV applied directional filter banks to suppress noise and new estimation kernels with a single-wave assumption [21]. However, the use of the Laplace operator may over-enhance the object boundary and the assumption of single-wave may be invalid when there are strong reflections. Elastography Software Pipeline (ESP) [22] used a series of filters for denoising and Gabor wavelets for inversion. ESP could recover refined anatomical details but required many parameter-tuning steps.

Non-Linear Inversion (NLI) algorithms based on Finite Element (FE) can provide good inversion [9, 23, 24]. However, the computation cost was high and the boundary condition settings could greatly influence the results [25]. Notably, by implementing a parallelized, subzone-based domain decomposition approach, the NLI algorithm proposed by the Dartmouth group avoids undue bias from the applied boundary conditions and high computational cost [26]. The method has been used to establish benchmark values for viscoelastic properties of the human brain [27]. Recently, data-driven algorithms using Deep Learning (DL) were proposed using the training set from FE simulation [28, 29]. These methods also relied on specific FE models for specific application scenarios, hindering their generalization performance. Traveling Wave Expansion (TWE) was first introduced to solve the inversion problem by using large fitting windows [30]. However, the inverse of the ill-conditioned matrix limited the spatial resolution that can be achieved. Therefore, a noise-resistant method with low-computational-cost, clearly delineated anatomical features and strong generalizability is required.

### Paper Contribution

In this study, we propose an optimization-based phase unwrapping method with Dual Data Consistency (Dual-DC) to simultaneously perform phase unwrapping and FFT. For the inversion of complex shear modulus, a TWE-based Neural Network (TWENN) algorithm is proposed. Detailed Simulation, phantom, and human studies were performed to demonstrate the pipeline's accuracy and reliability.

### Paper Organization

The remaining sections are organized as follows: In Section II, the proposed pipeline including Dual-DC and TWENN is presented. In Section III, the test datasets, evaluation metrics, implementation details of the proposed pipeline and comparative algorithms are introduced. Section IV presents detailed results followed by Section V that provides a thorough analysis of the results.

## II Methods

The proposed pipeline for processing MRE images includes an optimization-based phase unwrapping method with Dual-DC, and a Traveling Wave Expansion-based Neural Network (TWENN) modulus estimation to estimate complex shear modulus (Fig. 1).

### Optimization-based Phase Unwrapping with Dual Data Consistency

The displacement extraction of phase images acquired from MRE generally consists of two steps: phase unwrapping and computation of the principal components using FFT [7]. Here, an optimization-based method was proposed to perform these two tasks simultaneously.

#### Ii-A1 Problem Formulation

If \(\mathbf{U}^{*}=\mathbf{U}^{\prime}+\mathbf{i}\cdot\mathbf{U}^{\prime\prime}\) is the complex displacement field rotating at a frequency \(\omega\), representing the oscillating motion. In the process of image acquisition, at a specific time point \(t_{j}(j=1,...,J)\) within the motion cycle where the phase offset is \(\varphi_{j}=\omega t_{j}\), the image phase recorded by motion encoding gradient is \(\mathbf{\phi}_{j}\):

\[\mathbf{\phi}_{j}=\mathrm{Re}\big{(}\mathbf{U}^{*}\cdot\mathrm{e}^{i\varphi_{j}}\big{)} =\mathbf{U}^{*}\cdot\mathrm{cos}\big{(}\varphi_{j}\big{)}-\mathbf{U}^{\prime\prime} \cdot\mathrm{sin}\big{(}\varphi_{j}\big{)}, \tag{1}\]

Therefore, the image \(\mathbf{I}_{j}\) acquired at \(t_{j}\) is

\[\mathbf{I}_{j}=|\mathbf{A}|\cdot\mathrm{e}^{i\varphi_{j}}\cdot\mathrm{e}^{i\mathbf{\phi}_{ j}}, \tag{2}\]

where \(|\mathbf{A}|\) is the magnitude of image and \(\mathbf{\Phi}\) is the background phase of \(\mathbf{I}_{j}\). The purpose of displacement extraction is to obtain \(\mathbf{U}^{*}\) from \(\mathbf{I}_{j}(j=1\...J)\), which includes two tasks: computation of principal component and phase unwrapping. Two objective functions are proposed to solve them, respectively.

#### Ii-A2 Objective Function for Principal Components

Cross differences among the acquired images are used to remove the background phase. For images acquired at two temporal points \(t_{p}\) and \(t_{q}\), \((p,q)\in\{(x,y)|x<y,1\leq x,y\leq J\}\):

\[\begin{split}\frac{\mathbf{I}_{p}}{\mathbf{I}_{q}}&=e^{i \left(\Phi_{p}-\Phi_{q}\right)}\\ &=e^{i\left(\Phi^{\prime}\cdot\left\{\mathrm{cos}\left(\varphi_{ p}\right)-\mathrm{cos}\left(\varphi_{q}\right)\right\}+\mathbf{U}^{\prime\prime}\cdot \left[-\mathrm{sin}\left(\varphi_{p}\right)+\mathrm{sin}\left(\varphi_{q} \right)\right]\right\}\end{split} \tag{3}\]

Define the constant coefficients \(E_{pq}\) and \(F_{pa}\),

\[\begin{split} E_{pq}&=\mathrm{cos}\big{(}\varphi_{p} \big{)}-\mathrm{cos}\big{(}\varphi_{q}\big{)}\\ F_{pq}&=-\mathrm{sin}\big{(}\varphi_{p}\big{)}+ \mathrm{sin}\big{(}\varphi_{q}\big{)}\end{split} \tag{4}\]

The relationship between unknown \(\mathbf{U}^{*}\) and known \(\mathbf{I}_{j}\) is established so that the principal components are obtained. Thus, the first data consistency (DC) term of phase can be written as:

\[\mathcal{O}_{DC1}=\sum_{p,q}\left\|\frac{\mathbf{I}_{p}}{\mathbf{I}_{q}}-e^{i\left(\mathbf{ U}^{\prime\prime}\cdot E_{pq}+\mathbf{U}^{\prime\prime}\cdot F_{pq}\right)}\right\|_{2} \tag{5}\]
#### 2.1.3 Objective Function for Unwrapping

With Eq. (5) only, \(\mathbf{U}^{*}\) may still be wrapped. The absolute value of the wrapped phase gradient is much higher than that after unwrapping. Using the chain rule, the gradient of phase after unwrapping (true phase gradient) can be obtained directly from the wrapped phase [7]. Thus, data consistency is achieved between the true phase gradient and the gradient of \(\mathbf{U}^{*}\).

If the background phase is ignored, \(\frac{t_{j}}{|t_{j}|}\approx e^{i\mathbf{\phi}_{j}}\), by using the chain rule, the phase gradient is calculated from the wrapped phase,

\[\frac{\partial e^{i\mathbf{\phi}_{j}}}{\partial x}=i\cdot e^{i\mathbf{\phi}_{j}}\cdot \frac{\partial\left(\mathbf{U}^{\prime}\cdot\cos\bigl{(}\varphi_{j}\bigr{)}-\mathbf{U} ^{\prime\prime}\cdot\sin\bigl{(}\varphi_{j}\bigr{)}\right)}{\partial x} \tag{6}\]

\[\cos\bigl{(}\varphi_{j}\bigr{)}\frac{\partial\mathbf{u}^{\prime}}{\partial x}-\sin \bigl{(}\varphi_{j}\bigr{)}\frac{\partial\mathbf{u}^{\prime\prime}}{\partial x}=- \frac{i}{e^{i\mathbf{\phi}_{j}}}\frac{\partial e^{i\mathbf{\phi}_{j}}}{\partial x}.\]

In matrix form, we have

\[\begin{bmatrix}\cos\bigl{(}\varphi_{1}\bigr{)}&-\sin\bigl{(}\varphi_{1}\bigr{)} \\ \vdots&\vdots\\ \cos\bigl{(}\varphi_{j}\bigr{)}&-\sin\bigl{(}\varphi_{j}\bigr{)}\end{bmatrix} \begin{bmatrix}\frac{\partial\mathbf{u}^{\prime}}{\partial x}\\ \frac{\partial\mathbf{u}^{\prime\prime}}{\partial x}\\ \frac{\partial\mathbf{u}^{\prime\prime}}{\partial x}\end{bmatrix}\begin{bmatrix} \frac{\partial\mathbf{u}^{\prime}}{\partial x}\\ \frac{\partial\mathbf{u}^{\prime}}{\partial x}\end{bmatrix}\begin{bmatrix}\frac{ \partial\mathbf{u}^{\prime}}{\partial x}\\ \frac{\partial\mathbf{u}^{\prime}}{\partial x}\end{bmatrix} \tag{7}\]

\[\text{Let }\mathbf{H}=\begin{bmatrix}\cos\bigl{(}\varphi_{1}\bigr{)}&-\sin \bigl{(}\varphi_{1}\bigr{)}\\ \vdots&\vdots\\ \cos\bigl{(}\varphi_{j}\bigr{)}&-\sin\bigl{(}\varphi_{j}\bigr{)}\end{bmatrix} \begin{bmatrix}\frac{\partial\mathbf{u}^{\prime}}{\partial x}\\ \frac{\partial\mathbf{u}^{\prime}}{\partial x}\end{bmatrix}\begin{bmatrix}\frac{ \partial\mathbf{u}^{\prime}}{\partial x}\\ \frac{\partial\mathbf{u}^{\prime\prime}}{\partial x}\end{bmatrix} \tag{8}\]

Thus, the second data consistency term (DC) of phase gradient can be defined as

\[\begin{split}&\mathcal{O}_{DC2}=\left\lVert\frac{\partial\mathbf{U} ^{\prime}}{\partial x}-\mathbf{\mathbf{U}}^{\prime\prime}x\right\rVert_{2}+\left\lVert \frac{\partial\mathbf{U}^{\prime}}{\partial y}-\mathbf{\mathbf{U}}^{\prime\prime}y\right\rVert _{2}\\ &+\left\lVert\frac{\partial\mathbf{U}^{\prime\prime}}{\partial x}-\mathbf{ \mathbf{U}}^{\prime\prime}x\right\rVert_{2}+\left\lVert\frac{\partial\mathbf{U}^{ \prime\prime}}{\partial y}-\mathbf{\mathbf{U}}^{\prime\prime}y\right\rVert_{2}.\end{split} \tag{9}\]

#### 2.1.4 Optimization Using Dual-DC

The displacement extraction problem in MRE including FFT and unwrapping can be treated as an optimization problem (10) by combining (5) and (9) as an objective function. Here, the adaptive Momentum (ADAM) algorithm was used to solve this optimization problem.

\[\begin{split}\min\limits_{\mathbf{U}}\left\{\sum_{p,q}\left\lVert \frac{t_{p}}{t_{q}}-e^{i\left(\mathbf{u}^{\prime}\cdot E_{pq}+\mathbf{U}^{\prime\prime }\cdot E_{pq}\right)}\right\rVert_{2}+\right.\\ \left.\lambda\left(\begin{array}{c}\left\lVert\frac{\partial\mathbf{u}^{ \prime}}{\partial x}-\mathbf{\mathbf{U}}^{\prime}x\right\rVert_{2}+\left\lVert\frac{ \partial\mathbf{u}^{\prime}}{\partial y}-\mathbf{\mathbf{U}}^{\prime}y\right\rVert_{2}\\ \left.+\left\lVert\frac{\partial\mathbf{u}^{\prime\prime}}{\partial x}-\mathbf{\mathbf{U}}^ {\prime\prime}x\right\rVert_{2}+\left\lVert\frac{\partial\mathbf{u}^{\prime \prime}}{\partial y}-\mathbf{\mathbf{U}}^{\prime\prime}y\right\rVert_{2}\right)\right\} \end{split}\right\}\end{split} \tag{10}\]

where \(\lambda\) is a weighting parameter.

**Input**: Measured complex MR images \(\mathbf{I}_{j}\), Phase offset \(\varphi_{j}\) for each image.

**Output**: Main displacement components \(\mathbf{U}^{*}=\mathbf{U}^{\prime}+i\cdot\mathbf{U}^{\prime\prime}\)

1. Initialize the \(\mathbf{U}^{*}\) to be updated and set weighting parameter \(\lambda\) and maximum number of iterations \(maxIter\). Set number of iterations \(L=0\).
2. Calculate cross phase differences \(\frac{t_{p}}{t_{q}}\). (3)
3. Calculate coefficients \(E_{pq}\) and \(F_{pq}\) using \(\varphi_{j}\). (4)
4. Calculate phase gradient \(\mathbf{\mathbf{U}}^{\prime}_{x}\), \(\mathbf{\mathbf{U}}^{\prime}_{y}\), \(\mathbf{\mathbf{U}}^{\prime\prime}x\), \(\mathbf{\mathbf{U}}^{\prime\prime}y\) from MR images \(\mathbf{I}_{j}\). (8)
5. **repeat**
6. Update \(\mathcal{O}_{DC1}\). (5)
7. Update the gradient of \(\mathbf{U}^{*}\), \(\frac{\partial\mathbf{u}^{\prime}}{\partial x}\), \(\frac{\partial\mathbf{u}^{\prime}}{\partial y}\), \(\frac{\partial\mathbf{u}^{\prime\prime}}{\partial x}\), \(\frac{\partial\mathbf{u}^{\prime\prime}}{\partial y}\).
8. Update \(\mathbf{\mathbf{O}}_{DC2}\)(9)
9. Update \(\mathbf{U}^{*}\) using ADAM algorithm. (10)
10. \(L=L+1\)
11. **until \(L\geq maxIter\)**
12. **return \(\mathbf{U}^{*}=\mathbf{U}^{*}+i\cdot\mathbf{U}^{\prime\prime}\)**

### _TWEEN: Traveling Wave Expansion-based Neural Network Modulus Estimation_

The TWE model was used to generate noisy training data. It was also used to construct and train a complex value covariance neural network to mine the mapping from wavefield to wavenumber, and finally use multi-frequency multi-directional fusion to obtain the complex shear modulus.

Fig. 1: A flow chart of the pipeline of MRE image processing for displacement extraction, training data generation, structure of network and multi-data fusion. Complex displacement field was extracted from phase images using an optimization of an objective function with dual-DC. Training data was generated by TWE model for network training. The normalized complex wavenumber was estimated by TWENN. Complex shear modulus was obtained by combining multi-frequency and multi-direction data.
#### Ii-B1 Traveling Wave Expansion

For homogeneous, incompressible, isotropic, and viscoelastic material, at spatial location \(r\), the complex value displacement is \(u(\mathbf{r})\). It is expressed as the superposition of \(M\) traveling waves [30]:

\[\begin{split}& u(\mathbf{r})=\sum\nolimits_{m=1}^{M}w_{m}\\ &=\sum\nolimits_{m=1}^{M}a_{m}\cdot e^{i\cdot\left(k^{\prime}+i \cdot k^{\prime\prime}\right)\cdot\widehat{a_{m}}\cdot\mathbf{r}}\end{split} \tag{11}\]

where \(w_{m}\) is the \(m\)th traveling waves, \(a_{m}\) is the complex value amplitude of the displacement, \(m\) is the number of traveling waves, \(\widehat{\mathbf{n_{m}}}\) is the unit vector of the propagation direction of the \(m\)th traveling wave, \(k^{*}=k^{\prime}+i\cdot k^{\prime\prime}\) is the local complex wavenumber, \(k^{\prime}\) is the real wavenumber related to material elasticity, \(k^{\prime\prime}\) is imaginary wavenumber related to material viscosity. The complex wavenumber can be normalized by vibration frequency \(\omega\):

\[\begin{split}\overline{k^{*}}&=\frac{k^{*}}{\omega }=\frac{k^{\prime}}{\omega}+i\cdot\frac{k^{\prime\prime}}{\omega}\\ &=\overline{k^{\prime}}+i\cdot\overline{k^{\prime\prime}}.\end{split} \tag{12}\]

Thus, the aim of the inversion is to find the operator \(\mathcal{F}_{\theta}\colon u\to\overline{k^{*}}\).

#### Ii-B2 Training Data Generation

To construct a neural network capable of realizing the mapping operator \(\mathcal{F}_{\theta}\), training data set was prepared using TWE model. The training set was created by varying different numbers of traveling waves \(M\), propagation direction \(\widehat{\mathbf{m_{m}}}\), complex amplitude \(a_{m}\) and complex wavenumber \(k\). In this study, an isotropic patch of 3\(\times\)3 (2D) or 3\(\times\)3\(\times\)3 (3D) was used for the inversion of \(k\) at each location. The relatively small patch used is to increase the spatial resolution. To obtain \(\mathcal{F}_{\theta}\) with better noise robustness, noise \(\zeta(\mathbf{r})\) was added to the training set. Suppose \(\mathbb{C}(\mathbf{r})\) is the normalized, complex-valued Gaussian noise, \(b\) is the intensity of the noise, \(\zeta(\mathbf{r})=b\cdot\mathbb{C}(\mathbf{r})\).The detail of dataset generation rules is shown in Appendix. The training data \(d(r)\) at location \(r\) can be written as:

\[\begin{split}& d(\mathbf{r})=u(\mathbf{r})+\zeta(\mathbf{r})\\ &=\sum\nolimits_{m=1}^{M}a_{m}\cdot e^{i\cdot\left(k^{\prime}+i \cdot k^{\prime\prime}\right)\cdot\widehat{a_{m}}\cdot\mathbf{r}}+b\cdot\mathbb{C }(\mathbf{r})\end{split} \tag{13}\]

For 2D cases, the training data \(\mathbf{d}\) can be written as a 2D tensor as follows:

\[\mathbf{d}=\begin{bmatrix}d\left(r_{-1,-1}\right)&d\left(r_{0,-1}\right)&d \left(r_{1,-1}\right)\\ d\left(r_{-1,0}\right)&d\left(r_{0,0}\right)&d\left(r_{1,0}\right)\\ d\left(r_{-1,1}\right)&d\left(r_{0,1}\right)&d\left(r_{1,1}\right) \end{bmatrix} \tag{14}\]

Where \(r_{x,y}\) is a vector pointing from the center of the patch to the target pixel position \((x,y)\). Likewise, for 3D cases, \(\mathbf{d}\) will be a three-dimensional tensor.

#### Ii-B3 Covariance Complex Neural Network

Studies have shown complex value neural networks are more suitable for solving complex-valued problems [31]. Therefore, a complex-valued neural network was constructed to process these complex-valued wavefields. In a cascade structure, TWENNs had 9-layer fully connected complex-value neural networks, with the number of neurons in each layer being 60,50,40,30,20,15,10,6,1, respectively. No extra layers were used.

If the weight matrix of the linear layer is \(\mathrm{W}=\mathrm{A}+i\mathrm{B}\), and the complex input is \(\mathrm{h}=\alpha+i\mathrm{\beta}\), the output is \(\mathrm{Wh}=\alpha\alpha-\mathrm{B}\beta+i\left(\mathrm{A}\beta+\mathrm{B} \alpha\right).\) To preserve the phase angle, and to keep the magnitude within a certain range, a new complex activation function named _modSigmoid_ is designed [32]:

\[\mathrm{modSigmoid}(\mathrm{z})=\mathrm{Sigmoid}([\mathrm{z}]+\mathrm{a}) \mathrm{e}^{i\theta_{\mathrm{z}}}, \tag{15}\]

where \(\mathrm{z}\) is the complex input data, \(\mathrm{a}\) is a real bias to be learned, \(\theta_{\mathrm{z}}\) is the angle of \(\mathrm{z}\) and \(\mathrm{Sigmoid}(\mathrm{x})=\frac{1}{1+e^{-\mathrm{x}^{2}}}\).

Because the covariance can expose the information within the wave field better [33], a covariance term \(\mathbf{D}=\mathrm{vec}(\mathbf{d})\mathrm{vec}(\mathbf{d})^{H}\) is used as the model input, where \(\mathrm{vec}(\cdot)\) is an operator transforms a tensor to a column vector. \(\vec{k^{\prime}}\) and \(\vec{k^{\prime\prime}}\) are estimated separately by two parallel multi-layer full connection networks (Fig. 1). Two separate networks are used to estimate \(\vec{k^{\prime}}\) and \(\vec{k^{\prime\prime}}\). For each network, mean square error is used as a loss function for optimization. In this study, the network was trained using the ADAM algorithm. \(\mathcal{F}_{\theta}(\omega)\) was trained and obtained at different frequencies \(\omega\).

#### Ii-B4 Multi-frequency and Multi-direction Fusion

For the wavefield \(\mathbf{U}_{\omega m,d_{n}}^{*}\) obtained from MRE at different frequencies \(\omega_{m}(m=1,...,M)\) and in different directions \(d_{n}(n=1,...,N)\), the corresponding normalized complex wavenumber \(\overline{k^{*}}(\omega_{m},d_{n})\) is obtained using the trained network \(\mathcal{F}_{\theta}(\omega_{m})\).

For complex shear modulus \(G^{*}=G^{\prime}+iG^{\prime\prime}\), using corresponding principle [17],

\[G^{*}=\rho\frac{\omega^{2}}{(k^{\prime}-k^{*\prime\prime})^{2}}=\frac{\rho}{( \vec{k^{\prime}}-k^{*\prime\prime})^{2}}. \tag{16}\]

Thus, the storage modulus \(G^{*}\) and loss modulus \(G^{\prime\prime}\) can be estimated from the average of \(\overline{k^{*}}(\omega_{m},d_{n})\):

\[G^{\prime}+iG^{\prime\prime}=\frac{\rho}{\left(\frac{\sum_{m,n}\vec{k^{\prime }}(\omega_{m},d_{n})}{MN}\right)^{2}}. \tag{17}\]

This method does not do any pre-filtering to filter the noise, and the final modulus is filtered by a 3\(\times\)3 median filter.

**Input**: Wavefields of multi-frequency and multi-direction \(\mathbf{U}_{\omega_{m}d_{m}}^{*}\), information of spatial resolution and vibration frequencies \(\omega_{m}(m=1,...,M)\)

**Output**: Complex shear modulus \(G^{*}=G^{\prime}+iG^{\prime\prime}\)

1. Set the patch size, and intensity of the noise \(b\).
2. Generate training data based on spatial resolution and vibration frequencies using the traveling wave expansion model. (13)
3. Construct covariance complex neural network \(\mathcal{F}_{\theta}\). (14,15)
4. Training \(\mathcal{F}_{\theta}(\omega)\) using ADAM.
5. Calculate \(\overline{k^{*}}(\omega_{n},d_{m})\) using \(\mathcal{F}_{\theta}\).
6. Estimate \(G^{*}=G^{\prime}+iG^{\prime\prime}\) using multi-frequency multi-directional data fusion. (17)
7. **return \(G^{*}=G^{*}+iG^{\prime\prime}\)**

## III Validation and Verification

Wrapped phase images from simulation liver dataset and in
vivo brain MRE experiment were used for the validation of the proposed phase unwrapping method with Dual-DC. Modulus estimation with TWENN was validated using simulated phantom, brain, and liver wavefields. The proposed pipeline combining Dual-DC and TWENN was validated using phantom, brain, or liver MRE images. A summary of the data sets and methods used is shown in TABLE I. Results are also compared with those from state-of-the-art methods.

### Validation of Dual-DC phase unwrapping

#### 3.1.1 Simulated Liver Dataset

The artificially wrapped phase was generated using 2D wavefield data for publicly available simulated liver images ([https://bioqic-apps.charite.de/](https://bioqic-apps.charite.de/)). The displacements were normalized with a maximum phase of 4\(\pi\). Phases at 4 different temporal points within one cycle were generated. Complex Gaussian noise with different intensities was added to the complex image to evaluate the robustness of the algorithm.

#### 3.1.2 In vivo 3D Brain Dataset

The in vivo brain MRE images were acquired using a 3T scanner (uMR790, United Imaging Healthcare, Shanghai, China) with TR/TE=4000/65ms, MEG=40mT/m, resolution=3 mm \(\times\) 3 mm \(\times\) 3 mm, and 8 phase offsets [34]. The study protocol was reviewed and approved by the Institution Review Board of Shanghai Jiao Tong University.

#### 3.2.3 Algorithm Comparison

The performance of phase unwrapping method is compared with that from SR algorithm [12] and LBE [13, 35] that are commonly used in postprocessing MRE images [8]. Simulation results are compared with the Ground Truth (GT) while mean errors are calculated for validation. Interlayer continuity of the unwrapped phase is compared for brain MRE images.

### Validation of Modulus Estimation using TWENN

#### 3.2.1 Test Data Set

The test data set contained 10\({}^{7}\) examples for model evaluation. Wavefields with different signal-to-noise ratios (SNRs) and complex wavenumbers were generated using the TWE model, with a patch size of 3\(\times\)3 and a resolution of 3 mm \(\times\) 3 mm at 60Hz. SNRs were distributed uniformly from 12 dB to 38 dB, real normalized wavenumbers were uniformly distributed from 0.35s/m to 1.35s/m, while imaginary normalized wavenumbers were uniformly distributed between 0s/m and 0.28s/m. Estimation of \(\vec{k}^{\prime}\) and \(\vec{k}^{\prime\prime}\) at different SNRs and complex

\begin{table}
\begin{tabular}{c c c c c c c} \hline \hline
**Application** & **Data Set** & **Details of data** & **Data** & **Method** & **Methods for** & **Evaluation** \\  & & **Source** & **Proposed** & **Comparison** & **Metrics** \\ \hline \multirow{3}{*}{Phase} & Simulated & 2mm\(\times\)2mm; & \multirow{3}{*}{BIOQIC} & \multirow{3}{*}{Dual-DC} & SG & Mean error \\  & Liver & 42Hz; RO & & & & \\ \cline{1-1} \cline{4-7}  & \multirow{2}{*}{Brain MRE} & 3mm\(\times\)3mm\(\times\)3mm; & 3T & \multirow{2}{*}{DLE} & Interlayer \\  & & 50Hz; RO & scanner & & & \\ \hline \multirow{6}{*}{Modulus estimation} & \multirow{2}{*}{Test data} & 3mm\(\times\)3mm; & TWE & \multirow{2}{*}{TWEEN} & DI (2D) & Mean error \\  & & 60Hz; SS & model & & & \\ \cline{1-1} \cline{4-7}  & Simulated & 1.5mm\(\times\)1.5mm; & \multirow{2}{*}{COMSOL} & \multirow{2}{*}{k-MDEV} & \multirow{2}{*}{RMSE} \\  & phantom & 60,80,100Hz; SS & & & & \\ \cline{1-1} \cline{4-7}  & \multirow{2}{*}{Simulated liver} & 1.5mm\(\times\)1.5mm; & \multirow{2}{*}{BIOQIC} & \multirow{2}{*}{TWEEN} & MDEV & \\  & & 24, 28, 32, 36, 40, 44, 48, 52, & & & & \\ \cline{1-1} \cline{4-7}  & \multirow{2}{*}{BIOQIC} & \multirow{2}{*}{TWEEN (3D)} & MDEV & \multirow{2}{*}{RMSE} \\  & & 56, 60Hz; RO, PE, SS & & & & \\ \cline{1-1} \cline{4-7}  & \multirow{2}{*}{BIOQIC} & 2mm\(\times\)2mm; & \multirow{2}{*}{BIOQIC} & \multirow{2}{*}{LBE} & Match with \\  & & 20, 25, 30, 35, 40, 45Hz; & & & & \\ \cline{1-1} \cline{4-7}  & \multirow{2}{*}{Normal brain} & 200, PE, SS & & & & \\ \cline{1-1} \cline{4-7}  & \multirow{2}{*}{Normal liver} & & & & & \\ \cline{1-1} \cline{4-7}  & \multirow{2}{*}{Hepatic siderosis} & 2.7mm\(\times\)2.7mm; & \multirow{2}{*}{1.5T} & \multirow{2}{*}{MDEV} & \multirow{2}{*}{Regional mean} \\  & & 30, 40, 50, 60Hz; & scanner & & & \\ \cline{1-1} \cline{4-7}  & \multirow{2}{*}{Liver tumor} & & & & & \\ \cline{1-1} \cline{4-7}  & \multirow{2}{*}{Brain tumor} & 3mm\(\times\)3mm\(\times\)3mm; & 3T & \multirow{2}{*}{Dial-DC +} & \multirow{2}{*}{TWEEN} & \\  & & 30, 40, 50Hz; RO, PE, SS & scanner & & & \\ \cline{1-1} \cline{4-7}  & & & & & & \\ \hline \hline \end{tabular}
\end{table} TABLE I: A summary of data sets and methods used for comparison. 2D: 2D modelling inversion; 3D: 3D MOOLUS inversion; BIOQIC: DATASTS OR MRE processing pipeline were from [https://bioqic-apps.charite.de](https://bioqic-apps.charite.de); COMSOL: DATA was SIMulated using COMSOL (COMSOL AB, Stockholm, Sweden); 1.5T/3T scanner: MRE images were acquired from a 1.5T/3T scanner. Image RESOLUTION, Vibration FREQUENCY, AND THE MOTION ENCODING DIRECTION (RO: Read OUT, PE: Phase Encoding, SS: Slice Selection); ARE provolved IN THE DETIALS OF DATA.
wavenumbers were compared with those from conventional DI (fitting window=3\(\times\)3).

#### Iii-B2 Simulated Dataset

The size of the simulated phantom using COMSOL (COMSOL, Stockholm, Sweden) had two circular inclusions of 10 mm radius with complex shear moduli of 4 + i0.48 and 6 + i0.84 kPa. The complex shear modulus of background was 2 + i0.2 kPa. These values served as GT. A 28dB SNR of Gaussian noise was added to the wave images.

Publicly available MRE data for simulations of the liver and brain ([https://biooic-apps.charite.de/](https://biooic-apps.charite.de/)) were used for evaluation. Both were 3D, FE simulation models based on COMSOL. Model geometries were constructed from segmented images of a healthy human. Linear viscoelastic models were used and the tissue properties were assigned based on previously reported values [36, 37]. The data set contained a 3D multi-frequency and multi-direction wavefield with known \(G^{\prime}\) (GT). The liver dataset had a resolution of 2 mm \(\times\) 2 mm \(\times\) 2 mm with 4 frequency components (30, 36, 42, 48 Hz). The brain data set had a resolution of 1.5 mm \(\times\) 1.5 mm \(\times\) 1.5 mm, and 10 frequency components (24, 28, 32, 36, 40, 44, 48, 52, 56, 60 Hz).

Root mean squared error (RMSE) of estimated modulus \(G_{est}\)\({}^{\prime}\) and the structures were compared with k-MDEV and MDEV [38]. RMSE is defined as:

\[\text{RMSE}=\sqrt{\frac{\|\overline{G_{est}^{\prime}}-G_{GT}^{\prime}\|_{2}^{ 2}}{N}} \tag{18}\]

### Validation for the Proposed Pipeline

#### Iii-C1 Comparison with Other Pipelines

Both phantom and human MRE experiment images were used to evaluate the performance of the proposed pipeline for MRE image processing.

The proposed pipeline contains a cascading of phase unwrapping method using Dual-DC and TWENN for modulus estimation. Comparisons were made with two state-of-the-art multi-frequency MRE processing pipelines available on a public platform [38]: LBE and k-MDEV, FG and MDEV.

#### Iii-C2 Phantom MRE

For phantom validation, the multi-frequency phantom dataset containing four cylindrical inclusions was used [20]. The moduli estimated by each method were averaged in the z-direction. Regional mean values and ground truth of \(G^{\prime}\) and \(G^{\prime\prime}\) were compared.

#### Iii-C3 Brain MRE

Brain MRE images of a patient with a meningioma were acquired using a 3T scanner (uMR790, United Imaging Healthcare, Shanghai, China) using an electromagnetic actuator [34]. MRE was implemented using a single-shot spin-echo echo-planar imaging sequence (TR/TE=4000/65ms) with a motion-encoding gradient of 40 mT/m and 8 phase offsets. Contrast-to-noise ratio (CNR) of \(G^{\prime}\) with respect to the tumor region was used for evaluation:

\[\text{CNR}=\frac{2\left(\overline{G_{bkg}^{\prime}}-\overline{G_{tumor}^{ \prime}}\right)^{2}}{\sigma_{bkg}^{2}+\sigma_{tumor}^{2}}, \tag{19}\]

where \(\overline{G_{bkg}^{\prime}}\) and \(\overline{G_{tumor}^{\prime}}\) are the mean value of shear modulus, \(\sigma_{bkg}\) and \(\sigma_{tumor}\) are the standard deviation of shear modulus, for background tissue and tumor tissue, respectively. The tumor region was delineated from the T1 structural image, and the background neighborhood was delineated as the region of one tumor radius outside of the tumor.

A public data set of normal brain MRE images ([https://bioiqic-apps.charite.de/](https://bioiqic-apps.charite.de/)) was also used for evaluation by comparing the estimated modulus map and the structural image.

#### Iii-C4 Liver MRE

Liver MRE images from a healthy volunteer, a patient with a liver tumor, and a patient with hepatic siderosis & cirrhosis were acquired using a 1.5T pneumatic MRE system [39] (Magnetom Aera, Siemens, Erlangen, Germany, BIOQIC). CNR of \(G^{\prime}\) was used to evaluate the modulus estimation of liver tumor data. For normal liver data and iron-deposited cirrhotic data, the mean values of \(G^{\prime}\) of circled liver tissue were compared.

#### Iii-D4 Implementation Details for Dual-DC and TWENN

For all datasets unwrapped in this paper, the same hyperparameters of Dual-DC were used (learning rate = 0.005, \(\lambda=1000\), the number of iterations = 4000). The networks were trained using the ADAM algorithm with a batch size of 500 and a step number of 12000. Before the complex wave field was introduced into the network, complex normalization was performed based on the maximum absolute value of all data points. In this way, the absolute value of the input data was normalized between 0 to 1. Due to the fact that physiological noise such as vascular pulsation may pollute _in vivo_ measurement of brain and liver [40, 41], a relatively larger noise intensity \(b\) of 0.3 was added to the training sets, compared with that of 0.001 for the remaining non-living measurement

Fig. 2: (a) A comparison of the performance of the three unwrapping methods at varying noise levels. (b) The differences between the ground truth and the outputs of three methods. (c) Mean error at different noise levels.
dataset (all simulation data and phantom data). If the wavefield was continuous between imaging slices, a 3x3x3 patch was used. Otherwise, 3x3 patch was used.

The proposed methods were implemented with PyTorch 1.12.0 and CUDA 11.6 on an Ubuntu 20.04 LTS (64-bit) operating system equipped with an AMD Ryzen 9 5950x central processing unit (CPU) and NVIDIA RTX 3080Ti graphics processing unit (GPU, 12 GB memory).

## IV Results

### _Phase unwrapping test_

#### Iv-A1 Simulation test

Results of unwrapping simulation data showed the proposed optimization-based algorithm using dual-DC could unwrap phases with a complex Gaussian noise level up to \(\sigma=0.4\) without observing the offset of the background phase (Fig. 2a). However, SG did not perform unwrap at some complex edges without noise, and a large number of unwrapping failures were observed with increased noise levels (Fig. 2a, b). Although LBE had no obvious unwrapping failures, it introduced background offset even without noise (Fig. 2b). In contrast, an optimization-based algorithm using Dual-DC had better performance at different noise levels (Fig. 2c).

#### Iv-A2 In vivo dataset test

It was observed that SG was prone to unwrapping failure. The background offset introduced by the LBE could be seen in both coronal and sagittal views leading to an obvious 3D wavefield discontinuity. All these methods performed unwrapping at the transverse layers. However, an optimization-based algorithm

using Dual-DC still ensured continuity at coronal and sagittal views (Fig. 3). The absolute differences also showed that Dual-DC produced the smoothest phase.

### _Modulus estimation test_

#### Iv-B1 Simulation data set

A comparison of inversion between DI and TWENN with different complex wavenumbers and signal-to-noise ratios (SNR) in 2D wavefield (patch size=3x3) showed TWENN performed better than that of DI in estimating \(\vec{k^{\prime}}\) and \(\vec{k^{\prime\prime}}\). As SNR decreased, the advantages of TWENN over DI increased (Fig. 4a, b). Furthermore, the mean errors increased with the value of modulus, indicating an upper limit of estimation using TWENN. The upper limit, however, is higher than that of DI as shown in the figure.

Results of the estimated shear moduli showed TWENN could recover the inclusions better than other methods (Fig. 5). Values of RMSE showed TWENN was the best at estimating complex shear moduli of the inclusions in most cases (TABLE II).

#### Iv-B2 Simulation test of brain and liver

Based on the simulated complex wavefield of the brain and liver, it was observed that TWENN had the best estimation

Fig. 4: Distribution of mean errors of estimating \(\vec{k^{\prime}}\) (left) and \(\vec{k^{\prime}}\) (right) for different SNR values and normalized complex wavenumber.

Fig. 5: Comparisons of estimated (a) \(G^{\prime}\) and (b) \(G^{\prime\prime}\) maps using different methods based on phantom simulation. 28dB of noise was added to the simulated wave images.

Fig. 3: (a) comparison of displacement field extracted from brain MRE images. (b) The corresponding absolute differences of displacement at cross sections marked by the white dotted line.
accuracy with RMSEs of 3.34 and 2.06kPa (TABLE II). The other two algorithms did not perform as well as TWENN probably due to the strong reflection of waves or the use of derivatives, which might produce artifacts for boundaries (Fig. 6).

### _Performance Evaluation_

#### Iv-C1 Phantom MRE

Using multi-frequency MRE images of phantom, estimated \(G^{\prime}\) and \(G^{\prime\prime}\) values showed the proposed pipeline had the best performance in recovering the inclusion structures (Fig. 7a). Algorithms using the Laplace operator could not recover the structures well, probably because phantom images contained various noise and wave deflections. With longer wavelengths and lower SNR in regions with relatively high modulus, these algorithms were prone to underestimate the high modulus. The proposed pipeline produced the closest modulus estimation to the ground truth for both \(G^{\prime}\) and \(G^{\prime\prime}\) at most cases (Fig. 7b).

#### Iv-C2 Tumor MRE

For both brain and liver tumors, the proposed pipeline had better reconstruction of boundaries (Fig. 8a, b), with more details matching that of T1 images and the highest values of CNR (TABLE II). This showed the ability of the proposed pipeline to reconstruct the structural details.

#### Iv-C3 Normal brain MRE

For the normal human brain dataset, \(G^{\prime}\) values estimated from MDEV were relatively low (\(<\)1kPa). This was probably affected by noise. Both the k-MDEV and the proposed pipeline could estimate shear modulus maps which match the structure image. k-MDEV introduced significant boundary artifacts which is consistent with the results in Fig. 6. As indicated by the white arrow in Slice 1, the brain's midline was better reconstructed using the proposed pipeline than that of other algorithms (Fig. 8c). The white arrow of Slice 2 (Fig. 8c) showed k-MDEV estimated relatively low modulus (\(<\)0.3 kPa), but the proposed method provided a proper estimate (\(\sim\)1kPa). This dataset demonstrated the capability of the proposed pipeline to reconstruct anatomical features.

#### Iv-C4 Liver dataset test

For liver MRE (Fig. 8e, f), the estimated \(G^{\prime}\) values of the normal liver from k-MDEV, MDEV, and the proposed pipeline were 1.62, 1.31, and 1.97 kPa, respectively. These values were all in line with the reported normal range [42]. However, in the case of cirrhotic with iron deposition, the estimated \(G^{\prime}\) values were 1.70 and 1.09 for k-MDEV and MDEV, respectively. This was likely caused by the iron deposition that resulted in a low SNR. The value of 3.13 kPa estimated by the proposed pipeline fell within the reported range [42].

## V Discussion

A pipeline for estimating complex moduli from MRE images is proposed in this study. The processing pipeline contains two major steps: displacement extraction and modulus estimation. For displacement extraction, an optimization-based method with Dual-DC terms is used for both calculating principal components and phase unwrapping. For modulus estimation, a complex neural network using covariance as input is proposed based on the TWE model.

Since discontinuities of unwrapped phase images can introduce serious artifacts of the derived modulus distribution, accurate phase unwrapping is an important prerequisite. Complex or noisy wavefields are difficult to deal with using search-like or discrete greedy update phase unwrapping techniques, such as the SG algorithm. Laplacian-based estimation (LBE) can perform unwrapping well in noisy cases, but extra background offset noise could be introduced. The Dual-DC based continuous optimization can deal with unwrapping of complex wrapped and noisy phase without introducing extra phase offset. The proposed method for unwrapping uses phase gradient information, which is more adaptable to complex scenarios than discrete unwrapping. Compared to LBE, Dual-DC does not introduce new background noise. Phase images were unwrapped successfully for wavefield data from phantom, brain and liver, using the same set of hyperparameters. This verifies the generalizability of the proposed algorithm. In addition, the computation time for

Fig. 6: A comparison of estimated \(\bar{c}^{\prime}\) using different algorithms for simulated brain and liver. Regions were magnified to illustrate the anatomical features of \(\bar{c}^{\prime}\).

Fig. 7: (**a**) A comparison of estimated \(\bar{c}^{\prime}\) and \(\bar{G}^{\prime\prime}\) maps of gelatin phantom using different processing methods. (**b**) A comparison of mean values of \(\bar{c}^{\prime}\) and \(\bar{G}^{\prime\prime}\) of each method. The longitudinal axis is displayed in a logarithmic form.
MRE images of the whole brain (3 frequencies and 3 directions) was about 3 minutes.

For modulus inversion, where noise robustness, structural resilience and computational complexity impose contradicting constraints, TWENN manages to reach a tradeoff among these three metrics. In terms of noise resilience, TWENN does not use explicit differentiation computation, which also provides structural resolution benefits by reducing potential edge artifacts. Except for the LFE-type and NLI-type methods [9, 25, 23, 26], most of current algorithms rely on differential computation, which can be sensitive to noise. By adding noise to the training set, TWENN can simultaneously achieve both denoising and inversion, improving its performance at regions with low SNR. For example, TWENN can estimate the modulus of cirrhotic tissue at low SNR regions induced by iron deposition (Fig. 8f). TWENN has the potential to improve the clinical performance of liver MRE due to the fact that iron deposition is a significant cause of liver MRE failure [42].

Considering structural resilience, a small patch of 3\(\times\)3 or 3\(\times\)3 pixels is applied so that TWENN can still distinguish details of the image objects even if local homogeneity is assumed. Benefiting from the absence of differential calculations, no significant artifacts are observed at the structural edges. Compare with MDEV and k-MDEV, both of which apply differential operations, TWENN has better performance in preserving anatomical details on both simulated and in vivo liver and brain datasets (Fig. 6, Fig. 8).

If the estimation kernel is chosen in a trial-and-error fashion, it is difficult to ensure the noise robustness of small patches. Therefore, neural networks are used here to establish the nonlinear mapping between noisy wavefields and complex wavenumbers. Even training without noise, TWENN still has better noise robustness than conventional DI using a small kernel (Fig. 4). The dimension of the kernel also affects the inversion. The use of a small sized patch will set an upper limit to the estimated modulus. Using dejittering method, it has been shown that interslice phase inconsistencies introduced in signal acquisition can be removed, providing robust 2D inversions, especially in the abdomen. Here, with improved phase unwrapping, the smooth and physically accurate phase can result in better modulus estimation from 3D inversion.

TWENN uses neural networks to pre-learn wavefield structure information for inversion. The time of training data generation and network training using a desktop PC did not exceed 3 minutes. The trained network can be directly used to process MRE images with the same vibration frequency and resolution. In this study, the 3D multi-directional and multi-frequency inversion time of a whole brain did not exceed 15 seconds, showing its potential to be deployed clinically. This is because that all these methods use multi-frequency fusion, which is not used by other leading FE-based approaches like NLI. Furthermore, inverse FE-based algorithms require iterative FE computation updates [25, 37, 43] which may result in computations ten times or higher than TWENN, MEDV, and k-MEDV.

The modulus inversion of MRE is equivalent to estimating the wavenumber for an array of wave points, which can be considered as a variant of the Direction of Arrival (DOA)

Fig. 8: T1-weighted (T1w) images, magnitude images of MRE, \(G^{*}\) and \(G^{*}\) maps estimated from different algorithms were compared for (a) brain tumor, (b) liver tumor, and (c, d) two difference slices of normal brain images. Magnified views of the tumor region were also provided. Magnitude images of MRE, \(G^{*}\) and \(G^{*}\) maps estimated from (e) a healthy volunteer and (f) a patient with hepatic siderosis & cirrhosis were compared. Regions of interest were delineated with white line, and the corresponding values of mean \(\pm\) standard deviation of \(G^{*}\) were provided.
problem in array signal processing. Strategies commonly used in DOA solution problems [30] can be drawn on, where second-order moments can effectively reveal harmonic information. The possibility of estimating modulus from an array signal processing perspective was validated in our previous work [44].

From the perspective of generalizability, the training data set in TWENN is produced using a wave equation. Theoretically, any local wavefield can be represented using the TWE model, including reflected standing waves. Therefore, the training set can cover most of the wave propagation conditions, ensuring the generalizability of TWENN. Existing neural network training modulus estimation methods [28, 29] are usually trained using FE simulation for a limited number of scenarios, which could hamper their generalization performance. In addition, most of the elastography processing pipelines used various complex filters for either wave extraction or denoising [20, 21, 22]. However, the proposed pipeline does not use any filters other than a median filter with fixed parameters. Therefore, complicated parameter tuning procedures are prevented, improving its generalization performance.

Although multi-frequency and multi-directional information is used to improve inversion performance [2], TWENN also supports single-frequency and single-direction inversion. In addition to \(G^{\prime}\) and \(G^{\prime\prime}\) values, shear wave speed and penetrating rate that are closely related with stiffness and damping [2] can also be estimated.

Limitations of this study include assumptions of local homogeneity, isotropy, and linear viscoelasticity, which could be addressed by more realistic FE simulation based NLI [9, 23, 24, 25, 45, 46], and neural network based inversion [29, 47, 48]. Although TWENN can estimate a relatively smoother viscosity distribution from the noisy wavefield, only simple phantom cases were used to validate it. In this study, limited human imaging data were used for validation. Future work includes applying the proposed pipeline to a larger cohort of clinical data sets such as neurological and liver diseases.

The unwrapping method proposed in this paper is readily transferable to other fields where unwrapping is required. The inversion method can also be applied to other elastic imaging modalities, such as optical elastography and ultrasound elastography where wave equations apply. The TWE model can be modified in terms of specific anisotropic and dispersive material models, in order to solve more complex modulus inversion problems. Future work includes estimating properties of transversely isotropic and frequency-dependent materials within the TWENN framework. The potential of using this data-driven method for obtaining inverse operators to solve other inverse problems such as electrical resistance tomography needs to be further explored.

## Training data generation

The training data set was generated using the following parameter settings:

1. The number \(M\) of traveling waves 1,2,3,4,5,6,7,8 was set to be distributed at ratios of 6:4:3:2:1:1:1:1:1:1.
2. The complex amplitude \(|a_{m}|\) was uniformly distributed in [0,1] and angle(\(a_{m}\)) was uniformly distributed in [0,2\(\pi\)].
3. The unit vector of the wave propagation direction is \(\vec{n}_{m}^{*}=[x,y,z]\). 1. In 2D cases, \(x=\cos(\theta)\,,y=\sin(\theta)\,,z=0\), \(\theta\) was uniformly distributed in [0,2\(\pi\)]. 2. In 3D cases, \(z\) is uniformly distributed in [\(-1\),1], \(\theta\) was uniformly distributed in [0,2\(\pi\)], \(x=\sqrt{1-z^{2}}\cos(\theta)\), \(y=\sqrt{1-z^{2}}\sin(\theta)\).
4. \(k^{\prime}\) and \(k^{\prime\prime}\) were uniformly distributed within the preset range.

## Acknowledgment

We thank Prof. Ingolf Sack from Charite, Germany for providing part of the testing data sets and helpful discussions.
