Ultrasound elastography using a regularized modified error in constitutive equations (MECE) approach: a comprehensive phantom study

Many of the current techniques in transient elastography, such as shear wave elastography (SWE) assume a dominant planar shear wave propagating in an infinite medium. This underlying assumption, however, can be easily violated in real scenarios in vivo, leading to image artifacts and reconstruction errors. Other approaches that are not bound to planar shear wave assumption, such solutions based on the partial differential equation, can potentially overcome the shortcomings of the conventional SWE. The main objective of this paper is to demonstrate the advantages of the modified error in constitutive equations (MECE) formulation with total variation regularization (MECE + TV) over SWE in reconstructing the elastic moduli of different tissue-mimicking phantoms. Experiments were conducted on phantoms with inclusions of well-defined shapes to study the reconstruction of specific features relevant to practical applications. We compared the performances of MECE + TV and SWE in terms of quantitative metrics to estimate reconstruction accuracy, inclusion shape recovery, edge preservation and edge sharpness, inclusion size representation, and shear elasticity and contrast accuracies. The results indicate that the MECE + TV approach outperforms SWE based on several of these metrics. It is concluded that, with further development, the proposed method may offer elastography reconstructions that are superior to SWE in clinical applications.


Introduction
Ultrasound shear wave elastography (SWE) holds the potential for disease assessment and diagnosis and management such as staging liver fibrosis in the liver (Barr 2018) as well identification of malignant tumors, in breast (Athanasiou et al 2010, Cosgrove et al 2012, Gregory et al 2015, Denis et al 2015a, 2015b, and in thyroid (Gregory et al 2018(Gregory et al , 2020. Many of the current techniques in transient elastography assume a dominant planar shear wave propagating in an infinite medium, where the local elasticity values can be inferred from local shear wave velocities. In this scenario, the local shear wave velocity is determined from a mesh grid of time of flight (TOF) measurements (Nightingale et al 2003, Bercoff et al 2004, Song et al 2012. From the practical perspective, SWE methods are appealing as they can be implemented simply and computations can be carried out in near real-time. The main drawback of these methods is that their primary assumptions (e.g. propagating plane waves with known direction in an unbounded domain) can be easily violated in real scenarios in vivo. That assumption leads to unwanted artifacts in the reconstructed images (Nightingale et al 2003) or a biased estimation of modulus (Zhao et al 2011). The partial differential equation (PDE)-based methods relax the underlying assumptions of planar shear waves and unbounded mediums and, hence, can handle general waveforms, heterogeneous media, and general boundary conditions (Feissel and Allix 2007, Banerjee et al, 2013, Diaz et al 2015. Therefore, they work with a general elastic or viscoelastic initial-boundary-value problem.
The motivation for the present work is to develop a PDE-based solution to ultrasound elastography with the aim of overcoming the shortcomings of the conventional SWE. We recently applied a computationally effective iterative strategy for elastography called the 'modified error in constitutive equation (MECE)' (Ghosh et al 2017). In MECE the inverse problem is formulated as a nonlinear PDE constrained optimization problem with the objective consisting of energy and displacement functions (Banerjee et al, 2013, Feissel andAllix 2007). The MECE approach offers attractive advantages over other PDE-constrained optimization approaches including robustness to initial guess and local minima as well as the ability to accommodate incomplete knowledge of boundary conditions. The MECE approach was successfully evaluated in elasticity imaging in (Banerjee et al, 2013, Warner et al 2014. In this paper, we provide a quantitative comparison of the MECE approach with and without augmentation by a TV regularization functional, which we refer to as 'MECE + TV' and 'MECE' respectively, and an autocorrelation-based SWE algorithm using experimental phantom data. In clinical imaging, the shape and borders of emerging abnormalities in the soft tissue background carry valuable diagnostic information. With this in mind, in addition to studying the accuracy of elasticity estimation, here we pay particular attention to quantitative analysis of the shape and borders. Conducting a comparative quantitative study of the elasticity imaging methods mentioned above in vivo would be infeasible because the ground truth is not known a priori, and thus differences between reconstructions could not be given an independent physical verification. Likewise, ex vivo experiments would lack an inclusion with a controlled geometry for investigating the shape and border reconstructions characteristics of the reconstructed elasticity image. Therefore, we chose to conduct this study on well-controlled phantoms with known geometries and known shear modulus values. Many anatomical structures of interest have narrow or pointed extensions, where such morphology has diagnostic value. For example, the spiculations of breast masses are hallmarks of malignancy (Franquet et al 1993). For this reason, we used phantoms consisting of a soft gel matrix with inclusions in the form of a narrow layer (strip) and sharply pointed triangles and studied the reconstructions of shear modulus in these phantoms. The narrowing width of the triangular inclusion near the sharp vertex of the triangle serves to demonstrate the limit to which such fine structures can be captured by the applied method. Furthermore, these geometrical features invalidate many of the assumptions made in SWE. The goal is to compare the performance of MECE + TV and SWE in the reconstruction of different shapes and edges and evaluate the accuracy in the estimation of the shear modulus values. To assess the edge recovery, several performance measures have been studied, including edge preservation index (EPI) (Sattar et al 1997), zero normalized cross-correlation (Z-NCC) (Haralick and Shapiro 1992), Dice coefficient similarity index (Dice 1945). Moreover, for layered phantom, the width of the inclusion layer was estimated and the estimation error was reported.
The remainder of the paper is organized as follows. We first describe the theoretical background of the MECE + TV method in sections 2. The SWE method is presented in section 3. In section 4.1, we provide details of the experiments conducted in this work. Then, the study proceeds with applying MECE + TV and SWE toward estimating the shear modulus distribution from the experimental measurements in section 4.2. We provide conclusions and future directions in section 5.

Theoretical background: the MECE and MECE + TV approaches
In this section, we provide a brief description of the MECE formulation (Banerjee et al, 2013, Warner et al 2014, Ghosh et al 2017. We focus our descriptions on the MECE formulation for unknown boundary conditions (Diaz et al 2015) because the formulation naturally accommodates measured data from SWE experiments, as demonstrated in (Banerjee et al, 2013, Ghosh et al 2017. We begin by describing our model of interest; we then proceed to explain how we use the MECE formulation to invert for the material parameters in the model.

Mathematical model
In this work, we assume we have a compressible, isotropic two-dimensional (2D) elastic body undergoing plane strain deformations. We further assume that the excitation on the body is time-harmonic. The conservation of momentum equations for this situation is where σ is the Cauchy stress tensor, ρ is the density, ω is the angular frequency, and u is the displacement field. The stress tensor is related to the shear and bulk modulus G, B through a constitutive equation built into the formulation of the inverse problem (section 2.2). In order to invert (1) for the material parameters, we need to solve a boundary value problem. We have, however, neglected to provide boundary conditions because we typically do not have this information from ultrasound measurements. The MECE formulation that we describe in section 2.2 is able to handle missing boundary conditions (Warner et al 2014).
In this work, we model the deformation of tissue-mimicking phantoms, which are designed to have the same acoustic and mechanical properties as soft tissue. These phantoms are therefore nearly incompressible, so we set their bulk modulus to that of water. The only unknown parameter in our model is the shear modulus, and this is what we solve for in the inverse problem described next.

Inverse problem
The inverse problem of interest is to reconstruct the shear modulus, given measurements of the displacement field u. This inverse problem can be solved using the regularized MECE formulation described below In the above equations, σ is the stress tensor, σ dev is the deviatoric part of the stress tensor, ε is the strain tensor, T is an observation tensor that maps computed displacements to measured quantities, p is the hydrostatic pressure, e u is the volumetric strain, and u m is the measured displacement field. Notice that equation (4) implies a (constitutive) relationship between deviatoric stress, deviatoric strain, and shear modulus. Also, this equation relates pressure to the bulk modulus and the volumetric strain. If a zero constitutive error can be achieved, we would recover the conventional elasticity constitutive equations.
The optimization problem described in equation (2) has an objective function with three components: an error in constitutive equation (ECE) term, a data misfit term, and a total variation regularization (TV) term. The ECE term is used to relax how strongly the data satisfies the constitutive equation. The regularization term is used to impose smoothness constraints on the shear modulus field. A TV regularization is chosen because it allows for sharp jumps in the reconstructed parameters while removing oscillations due to noise in the data (Vogel 2002).
The objective function contains various hyperparameters. The data misfit hyperparameter controls how closely the predicted displacement field matches the measurements . The data misfit term also contains a weighting tensor that can be used to place more weight on certain components of the displacement field. This is useful for ultrasound data, which has a much higher sampling rate in the axial direction than in the lateral direction. The regularization hyperparameter controls the smoothness and sharpness of the jumps in the reconstructed fields. The regularization term also contains a hyperparameter that is used to make sure that the derivative of the regularization term is well defined when ∇G = 0. Other quantities used in equation (4) are the strain displacement relationships, and the deviatoric and hydrostatic stresses and strains which are defined as p := 1 3 tr (σ) (10) Here, I. is the identity tensor. The optimization problem defined in equation (2) contains the function space U, which is an admissible space of displacements, S (u) which is an admissible space of stresses (i.e. symmetric tensors), and Z, which is an admissible space of shear modulus (i.e. positive and bounded). We solve this optimization problem by first differentiating the objective function with respect to the independent variables to obtain a system of nonlinear equations, which are discretized with the finite element method (FEM), and minimized with the alternating directions minimization algorithm (ADM) (Diaz et al 2015), and a hybrid algorithm which combines the ADM algorithm and a quasi-newton optimization algorithm. An in depth mathematical analysis of the MECE approach, which highlights some of its most salient features such as robustness to initial guess, treatment of missing boundary conditions, and resilience to resonances, is presented in (Diaz et al 2015)

Hyperparameter selection
The MECE objective function has three hyperparameters (β, κ, η) that need to be selected to solve the optimization problem defined in (2). We fix η = 10 −7 to make sure that the derivative of the TV regularization term is well defined when ∇G = 0. The other two hyperparameters β, κ are selected with variations of the error balance technique (Banerjee et al, 2013, Warner et al 2014 When we use just the ADM algorithm to minimize the unregularized MECE objective function, then we use the standard error balance technique defined in (Banerjee et al, 2013, Warner et al 2014 to select κ. When we use the hybrid algorithm to minimize the regularized MECE objective function, then we use the modified version of error balance technique to select κ, and β.

TOF SWE approach
Here, we briefly describe the SWE technique. For a shear wave propagating in an infinite medium, the shear wave velocity is proportional to the square root of the shear modulus. Local shear wave velocities are calculated using the TOF technique. In TOF, for each point in the imaging plane, a region of interest is selected and based on a known shear wave propagating direction, time delay, ∆t of wave trajectory is calculated to estimate the wave velocity c s . A bandpass filter with an effective bandwidth between 65 and 500 Hz is used to eliminate unwanted frequency components of the shear waves. Moreover, a Tukey time window technique with a parameter equal to 0.25 is applied to impose zero displacement at the beginning and the end of the displacement profile (Song et al 2012). The displacement profiles were interpolated by a factor of five before cross-correlation. The multi-directional filtering method reported in (Song et al 2014) provides a robust 2D shear wave map reconstruction for any arbitrary angle in the imaging plane.

Experiments description
Experiments were conducted on gel phantoms with agar inclusion. The inclusion material included agar material (Sigma-Aldrich, St. Louis, MO) with a concentration of 1.4% by volume in a mix of distilled water (50%) and evaporated milk (50%). Cellulose particles (Sigma-Aldrich, St. Louis, MO) with size 50 µm were also added with a concentration of 2% by volume to provide adequate ultrasonic scattering. Gelatin was used as the background material, using 300 Bloom gelatin and glycerol (Sigma-Aldrich, St. Louis, MO) with a concentration of 5% by volume. Potassium sorbate (Sigma-Aldrich, St. Louis, MO) at a concentration of 1% by volume was added as a preservative. All phantoms were fabricated in blocks of 80 × 80 mm gelatin-based material. Three-layered phantoms with a layer thickness of 3 mm, 5 mm and 8 mm were made. Moreover, one triangle inclusion phantom was made with a base of 11 mm and a height of 17 mm located at 20 mm from the phantom's top surface. The acoustic attenuation coefficient was 0.14 dB cm −1 -MHz for the background and 0.48 dB cm −1 -MHz for the inclusion materials. Mechanical testing was performed using a commercially available DMA device (RheoSpectris™ C500, QC, Canada) on three 4.5 cm long cylindrical tube samples (9 mm inner diameter) from each of the same gelatin and agar materials that were used for the phantom construction. The shear modulus of the phantom background and inclusion were determined to be 8.5 kPa and 17 kPa, respectively. These values were treated as the independent ground truth for the shear modulus of the phantom and were used later for comparison with our ultrasound measurements. A Verasonics V-1 programmable ultrasound machine (Verasonics Inc. Redmond, WA) equipped with a linear array transducer L7-4 (Philips Healthcare, Andover, MA) with a center frequency of 5 MHz was used for both phantom imaging and creation of the acoustic radiation force (ARF) push beams. For the ARF push, we used a tone-burst with a duration of 600 µs, consisting of 3 consecutive tone-bursts of 200 µs each. The push beam was followed by 3 angles (−2 • , 0 • , 2 • ) compounded plane-wave imaging. The imaging pulse repetition frequency (PRFd) was 3.3 kHz operating at 150 V peak to peak voltage; the length of data was 15 msec. Three experiments were conducted on each of the layered phantoms with different distances between the push location and layer inclusion. Four experiments were performed on the triangle phantom, two experiments per orientation with two different distances between the push location and the inclusion. In each experiment, a single ARF push was focused at a depth corresponding to the center of the inclusion; however, the push beam was applied at different lateral locations outside the inclusion. The sampling interval in the axial direction was set at λ/8, while the lateral sampling interval was set at λ, where λ is the wavelength corresponding to the center frequency of the probe (5 MHz). The displacement was then calculated for each point inside of the imaging plane using a 2D autocorrelation technique on the in-phase and quadrature (IQ) data (Loupas et al 1995). Eight axial samples were replaced by their average value to obtain a uniformly sampled velocity map in both axial and lateral directions. A schematic of the phantoms' geometry and location of pushes on different phantoms are shown in figure 1. It is noteworthy that in each experiment, we applied only one push beam at a time. The push beam location was changed between the experiments as shown in figure 1.

Methodology
Here, we present the reconstruction of the shear modulus using both MECE + TV and SWE methods and compare the results in the same region of interest. We use the registered geometry and location of inclusion within the phantoms on the B-mode image as the reference information to help compare various modulus maps (see figure 2).

MECE data pre-processing
Before processing the measured time-domain displacement field with the MECE formulation, we had to transform it into the frequency domain. This transformation was done with a fast Fourier transform (FFT) algorithm to obtain displacement fields at a broad range of frequencies. From these, nine displacement fields at frequencies between 104 and 520 Hz with 59 Hz increments, were selected and used for the inversion. This range is chosen by looking at the amplitude spectrum of the displacement field and selecting the frequencies that had a significant amount of energy. The lower bound of the frequency was chosen to be 104 Hz based on the fact that this frequency is larger than the inverse of the total time of the time-domain simulation.

Inversion with the MECE algorithm
We solve the inverse problem for the shear modulus using the axial displacement field because they are the only measured component. The components of the weighting tensor corresponding to the axial displacements are selected to be large relative to the components corresponding to the lateral displacements to account for the missing lateral measurements. We fix the density to 1000 kg m −3 , and the bulk modulus to that of water (2.2 GPa). We implemented the MECE inversion approach using our in-house finite element code written in C++ and used Linux workstations with multi-core Intel processors and 128 GB of RAM for all the solutions presented in this work.
The domain for the inversion is chosen to be a rectangular region within the entire field of view where the data was measured. Thus, the inversion domain boundaries are selected to exclude regions with poor measurements, such as the region containing the ARF push beam. This inversion domain is discretized with four-node bilinear quadrilateral finite elements. The mesh resolution is chosen so that there are at least 10 degrees of freedom per the wavelength corresponding to the highest frequency used for the inversion. The displacements are represented with piecewise bilinear shape functions, and the elastic modulus fields are represented with piecewise constant shape functions. The measured displacement field is interpolated onto the finite element mesh for the inversion.
The inverse problem is solved with and without TV. The unregularized MECE objective function is minimized with the ADM algorithm (Diaz et al 2015), while the regularized MECE objective function is minimized with the hybrid algorithm, which is a combination of the ADM and reduced MECE algorithms. For both algorithms, a constant initial guess of 5 kPa was used for the shear modulus. The ADM algorithm was terminated when the relative change in the MECE objective function drops below 10 −2 . For the hybrid algorithm, a loose tolerance of 10 −2 was used for the ADM algorithm, while the reduced MECE algorithm was terminated when any one of the following criteria is met: • The relative change in the objective function drops below 10 −4 .
• The norm of the gradient drops below 10 −5 .
• The norm of the change in the elastic modulus drops below 10 −10 .

Results and discussion
We evaluated the accuracy of the reconstructions by comparing the results to their known values with the following two criteria: (1) geometry of the reconstructed inclusion, and (2) the mean value of the shear modulus over the inclusion region and the background. Multiple geometrical parameters were considered. To assess the edge recovery several performance measures were studied, including edge preservation index (EPI) (Sattar et al 1997), ZNCC (Haralick and Shapiro 1992), and Dice coefficient similarity index (Dice 1945). Moreover, for the layer phantom, the width of the inclusion layer, as well as the associated normalized error, were estimated. Furthermore, the decay rate of modulus between inclusion and background is reported. Finally, to compare the elasticity values, the statistical evaluation of the shear modulus and contrast values are presented for different methods. Due to the noise in the MECE method (figures 3-7) it was excluded from the thresholding based metrics (i.e. Dice coefficient index, EPI, decay rate and width calculations) as such noise tended to saturate the thresholding described in section 4.3.2.2 and produce unreasonable binarized images.

Qualitative assessment
Using MECE, MECE + TV, and SWE methods, the reconstructed shear modulus maps for layered phantoms with the agar-based inclusion widths of 3 mm, 5 mm and 8 mm are presented in figures 3-5, respectively. The known boundaries have been identified using solid lines in the aforementioned figures.  The SWE reconstructed maps of modulus are illustrated in figures 3(g)-(i), 4(g)-(i) and 5(g)-(i) for different locations of ARF. These figures show that the general features, such as the position and geometry of the inclusions, are preserved to some extent. However, some anomalies, such as flattening of inclusion areas, are degrading in the SWE-based reconstructions.
Figures 3(j)-(l), 4(j)-(l) and 5(j)-(l) show the cross-sectional plots of the reconstructed shear modulus map of the phantoms with 3 mm, 5 mm and 8 mm layer inclusion along dashed line shown in figures 3(a)-(f), 4(a)-(f) and 5(a)-(f), for the MECE, MECE + TV, and SWE approaches, respectively. It is clear that although the field is fluctuating in the case of MECE, the inclusion is always stiffer than the background in all experiments using different methods. Moreover, the decay rate of shear modulus values from inclusion to the background region is slower in SWE compared to MECE + TV method, which indicates that the sharpness of the edges is better preserved by the latter method. The quantitative analysis of the shear modulus decay rate (edge sharpness) is given in subsection 4.3.3.1.
Many anatomical structures of interest have pointed extensions, where such extensions have a diagnostic value; for example, the spiculations of breast masses are a hallmark of malignancy (Franquet et al 1993). To assess the accuracy of the proposed method in the reconstruction of sharp and pointed edges, we study a phantom with triangular inclusions at two different orientations (see figures 6 and 7). The superimposed solid lines indicate the known boundaries of the inclusion.

Quantitative assessment
To evaluate the performance of shear modulus reconstruction, several quantitative assessment measures, including Zero-Normalized Cross-Correlation (Z-NCC), Contrast-to-Noise Ratio (CNR), Dice coefficient, modified EPI as well as reconstructed width of inclusion and decay rate of shear modulus for layer phantoms and the triangle phantom, are presented in this section.

Zero-normalized cross-correlation (Z-NCC)
Z-NCC is a known measure of template matching in computer vision (Haralick and Shapiro 1992). Here, it is used to compare the accuracy of reconstruction of the shear modulus map using SWE, MECE, and MECE + TV, with the known binary version of the inclusion (mask). This measure shows the correlation between the shear modulus reconstruction map and the ground truth. Therefore, a higher correlation value indicates a higher similarity between the reconstructed shear modulus image and the inclusion mask. It is defined as follows where I(m, n) is the reconstructed shear modulus image using MECE, MECE + TV or SWE. P is the phantom's binary inclusion mask, wherein pixels inside the inclusion mask are set to one and pixels outside the inclusion are set to zero.Ī m1,n1 is the average of I in the inclusion mask region,P is the mean value of the inclusion mask. m and n denote the pixel position and m 1 and n 1 the offset or lag of P relative to I. Similar to statistical correlation, values of γ(m 1 , n 1 ) that are closer to one indicate greater similarity between P and I, while values closer to zero indicate less similarity. As a complement to Z-NCC, we also calculated the CNR for each image. The CNR provides a measure of noise relative to the contrast provided by the inclusion. The CNR is defined as follows (Van Wijk and Thijssen 2002, Dumont et al 2016) Here, I inc and I back are the mean pixel values in the inclusion and background respectively and std(I inc ) and std(I back ) are the standard deviations of the pixels in the inclusion and background, respectively. Using an identical computational domain for MECE, MECE + TV and, SWE, the mean value of Z-NCC and CNR, averaged for different push beam locations, for each phantom are computed and plotted in figure 8. This figure shows that values closer to one are observed for MECE and MECE + TV compared to SWE for the triangular phantom. The difference between the mean of Z-NCC in MECE + TV and SWE are 0.31 and 0.26 for the triangle phantom at the two different orientations, which indicate that the results of reconstruction by MECE + TV are much closer to the inclusion mask than the reconstruction by SWE. In other words, MECE + TV outperforms SWE in the case of the triangular phantom. However, the Z-NCC values for the layer phantoms are almost the same for all reconstruction methods. This result is consistent with the fact that the layer appears wider than its original width in all reconstruction methods, hence the correlation with the ground truth in equation (12) results in large but relatively similar values.
Likewise, the CNR was also higher for MECE in SWE for the triangular inclusions. The mean differences were 2.26 and 2.21, respectively, for the two orientations. For the layer phantoms, the CNR for MECE + TV was nearly the same as SWE. Note that MECE by itself has a lower CNR than MECE + TV for all inclusions, consistent with the noise observed in it in figures 5-7.

Sørensen-Dice coefficients
The Sørensen-Dice (or Dice for short) coefficient (DSC) is a measure used for comparing the similarity of two images (Dice 1945). Here, the main reason to evaluate shear modulus reconstruction algorithms by the Dice coefficient is quantifying the performance of shape recovery in the phantoms, i.e. recovering the inclusion shape regardless of its value. Given two binary image of, (binarized image of the reconstructed shear modulus) and the binary image of the inclusion, (ground truth), DSC is defined as where |I b | and |P| are numbers of non-zero elements in each image and |I b ∩ P| is the number of none-zero elements in in which the corresponding element in is also none-zero. The ground truth image is binarized using simple thresholding. However, the binarizing process for I b is performed using three different levels of thresholds: 70%, 80% and 90% of the maximum amplitude of reconstructed shear modulus for each horizontal row of the reconstructed shear modulus image.
We computed and compared the Dice coefficients between shear modulus reconstruction and within the computational domain shown in figures 3-7 using MECE + TV and SWE methods. Figure 9 shows the mean of the Dice coefficient (averaged for different locations of the push beam and different values of thresholds) for each phantom. Although the overlap between Dice coefficients can be observed in figure 9 between MECE + TV and SWE, a paired t-test suggested the average difference between Dice coefficients to be statistically significant (p < 0.01; mean difference 0.15), supporting the observation that MECE + TV reconstruction method provides higher similarity to the ground truth than the SWE method. The difference between the mean of the Dice coefficients in MECE + TV and SWE methods are 0.11, 0.12, 0.17, 0.14, 0.24 for the 3 mm, 5 mm, 8 mm layered phantoms and for the first and second orientations of the triangle phantom, respectively.

Edge preservation of triangle phantom
To assess the performance of SWE and MECE + TV in edge preservation, we utilized an edge recovery measure that we call 'modified edge preservation index (MEPI)' . MEPI calculates the maximum value of zero-normalized cross-correlation for the binarized Laplacian-filtered versions of the inclusion ground truth and the reconstructed shear modulus images. MEPI can be considered as an extension of conventional EPI measure reported in (Sattar et al 1997). Figures 10(a)-(c) shows the recovered edge of the inclusion triangle and that of the reconstructed shear modulus map using MECE + TV and SWE. A visual inspection of this figure shows that the edges of the triangle inclusion using MECE + TV have been recovered more successfully compared to the corresponding edges obtained by SWE. The MEPI measure is defined as follows: where std(·) denotes the standard deviation of the reconstructed image andP E andĪ E are the mean pixel values of the respective image. A good match between the edge pixels of the reconstruction by MECE + TV or SWE and the ground truth results in a large MEPI. The MEPI is computed for the MECE + TV and SWE for the entire triangle inclusion and in the region of interest (indicated by a red square in figures 10(a)-(c)) containing the tip of the triangle. The reason for selecting this region of interest is that reconstruction of sharp and narrow structures are often challenging in elastography. Figure 10(d) compares the mean of MEPI, where the mean is calculated for different locations of the push beam (mean(MEPI)), in the triangle phantom at two orientations. The mean(MEPI) values of MECE + TV are consistently larger than those of SWE, particularly for orientation 1, demonstrate the superior performance of MECE + TV compare to SWE.

Estimation of inclusion width in layer phantoms
Another quantitative metric used to evaluate the quality of the reconstructions is the accuracy of the reconstruction method in preserving the geometry of the target. For this purpose, we computed an average width (mean(width)) of the inclusion in the layer phantoms as well as the average values of the normalized  error on the estimated widths (mean (NE)), defined as the average of the difference between the width of the layer in the reconstruction and that of the ground truth divided by the layer width in the ground truth. The averaging is performed on the data acquired for multiple push beam locations and multiple thresholding values (i.e. 70%, 80% and 90% of maximum amplitude in image). Figures 11(a) and (b) show comparisons between (mean(width)) and mean(NE) values using MECE + TV and SWE methods for different layer phantoms. The mean(NE) in MECE + TV were 18.62%, 9.52%, 11.15% for 3 mm, 5 mm and 8 mm layered phantoms, respectively; while the mean(NE) values using SWE were 36.98%, 18.28%, 14.65% for 3 mm, 5 mm, and 8 mm layer phantoms, respectively.

Edge sharpness
Another metric for comparison is the sharpness of the reconstruction at the inclusion-background interface of the layer phantoms. For this purpose, we evaluate the slope, or the decay rate, of the reconstruction at such interfaces. A higher decay rate implies a sharper, hence a more accurate, reconstruction of the interface. Figure 11(c) shows the mean of the decay rate of modulus at the interface between inclusion and background when the averaged is calculated for different locations of pushes in each layered phantom. It can be observed that the decay rate of the modulus is always larger in MECE + TV than in SWE, demonstrating that the MECE + TV reconstruction has preserved the sharpness of the edge at the inclusion-background interface more accurately.

Shear modulus and contrast values
The average shear modulus are computed for the inclusion and background of the phantoms within the computational domain shown in figures 3-7 using MECE, MECE + TV and SWE (See figures 12(a) and (b)). The average shear modulus for the background and inclusion were 8.9789 and 16.5977 kPa, respectively, using MECE, and 9.2473 and 16.07 kPa, respectively, using MECE + TV for all experiments. These values are in agreement with the known shear modulus values of the materials used in the experiments (8.5 kPa and 17 kPa) that were measured independently by the DMA device (RheoSpectris™ C500, QC, Canada) on samples of these materials. The observed variation in shear modulus is expected because of different levels of measurement and modelling errors in various experiments.
Another important parameter is the contrast. The Weber contrast is defined as the ratio of the difference between the inclusion and background average moduli to the average modulus of the background (figure 12(c)). The normalized error is defined as the difference between the estimated shear moduli in the reconstruction and the independently measured shear modulus of the inclusion and background materials. The normalized error of the shear modulus in the background and the inclusion are reported in figures 12(d) and (e), respectively. Further, figure 12(f) shows the normalized error of the contrast parameter, defined as the difference between the estimated contrast in the reconstruction and the true contrast (calculated using the independent elasticity measurements) divided by the true contrast.

Discussions and conclusions
In this paper, we compared two elasticity imaging methods: the traditional SWE and a newly developed full-wave inversion approach, MECE, and its regularized version, MECE + TV. The SWE method is based on the assumption of propagation of the plane shear wave in the object, whereas the MECE approach is not restricted to plane shear waves. In general, biological tissues may not support plane shear waves due to their complex geometry and heterogeneity. Therefore, it is expected that a full-wave reconstruction such as MECE would perform better than SWE. Based on this concept, the motivation for this paper was to explore and compare the performance of these two methods. Conducting a comparative study in vivo would be difficult because the ground truth is not known a priori. Therefore, we chose to conduct this study on well-controlled phantoms with known geometries and elasticity values. Two types of phantoms were designed for this purpose, a set of phantoms each with a layer of agar-based inclusion of different width in a uniform gel matrix, and a phantom with a triangular shape inclusion. Using these phantoms, we were able to test the performances of SWE, MECE, and MECE + TV based on the inclusion size, the sharpness of the inclusion edges, preservation of the geometry, and accuracy of shear modulus values and contrast.
In addition to qualitative evaluations, we used several quantitative measures to compare MECE, MECE + TV, and SWE. In particular, we provided metrics evaluating overall image quality and reconstruction accuracy (ZNCC, Contrast and Elasticity), inclusion shape recovery (Dice coefficient and layer width estimation) and inclusion edge recovery and sharpness (MEPI and Edge Decay Rate): Reconstruction accuracy: image quality and reconstruction accuracy were generally improved in MECE + TV compared to SWE, as demonstrated by the increased mean(ZNCC) and mean(CNR) results summarized in figure 8, showing that the reconstructions of the triangular phantom by MECE and MECE + TV have much higher correlation and CNR with the ground truth than the reconstruction using SWE. However, for the cases of layer phantoms, the MECE or MECE + TV do not show an improvement over SWE. This is expected because the SWE reconstruction of the layer inclusion generally results in a wider layer, for which in turn results in a higher correlation with the ground truth and less variation at the boundary of inclusion due to a smaller rising edge. That is, although the SWE reconstruction does not provide an accurate representation of the layer width, the cross-correlation and CNR are higher. Therefore, for the layer phantoms, a higher mean(Z-NCC) should not be interpreted as a more accurate reconstruction.
Shape recovery: shape recovery was evaluated using the Dice coefficient. The results shown in figure 9 showed that for all phantoms, the shape of the reconstructed inclusion is more similar to the known geometry of the phantom, i.e. the ground truth when the reconstruction is done by MECE + TV approach than SWE. It should be noted that since the Dice coefficient is based on the binary version of the image, the results do not apply to the elasticity values recovery.
Edge preservation: the mean of the modified edge preservation index (MEPI) was used to compare the preservation of the edges of the triangular inclusion in the MECE + TV and SWE reconstructions. Visual and qualitative inspection of the shape of the triangular inclusion shown in figures 10(a)-(c) demonstrates that MECE + TV outperforms SWE in preserving the triangle borders, particularly around the sharp corner at the bottom. A quantitative assessment of edge presented in figure 10(d), shows higher values of mean(MEPI) parameter for MECE + TV method than SWE for the region A around the bottom sharp corner of the triangular inclusion. We speculate that the failure of SWE in capturing the sharp corner is because the plane shear wave cannot be supported in areas with sharp transitions. In many clinical applications, borders of lesions carry valuable diagnostic value. For example, sharp extensions, or spiculations, of breast masses may indicate malignancy. The above results suggest that MECE + TV may be advantageous in such applications.
Size estimation: delineation of narrow structures in SWE elasticity reconstruction is often challenging. Figure 11(b) compares the normalized error in estimation of the width of the layer inclusions in MECE + TV and SWE reconstructions. Many anatomical structures, such as nerves and tendons, are narrow and long. The above results suggest that MECE + TV can improve the delineation of such organs in elasticity imaging.
Edge sharpness: another important factor in the delineation of the borders of narrow anatomical structures, such as nerves, blood vessels, and tendons, as well as delineation of the boundaries of isolated lesions is edge sharpness. The decay rate measure presented in figure 11(c) shows that the MECE + TV approach consistently preserves the sharpness of the layer inclusion boundaries with a higher decay rate than SWE.
Elasticity estimation and contrast: figures 12(d) and (e) showed that in most cases MECE reconstruction has a lower normalized error in elasticity value estimation than SWE. Another important factor is the elasticity contrast. In clinical practice, the contrast in elasticity is often used as a diagnostic biomarker. Our results in figure 12(f) show that MECE (and MECE + TV in most cases) are more accurate in estimating the elasticity contrast. These results suggest that MECE and MECE + TV offer higher diagnostic value than SWE.
There are some limitations to our study. We limited our study to only phantoms, mainly because the phantom material elasticity can be measured independently and reliably and used as a ground truth. However, phantom experiments such as these do not recreate several aspects relevant to in vivo experiments, such as low SNR and complicated elastic heterogeneity within the imaged domain. Additional studies are required to investigate how these results may translate to in vivo tissue measurements. The inclusions in our phantoms were uniform in the third (out of the page) dimension, thus effectively act as 2D inclusions. Studying phantoms with inclusions that are limited in all dimensions may be more realistic and provide information about the effect of inclusion volume on the elasticity reconstruction. Phantoms with spherical inclusions of varied diameters, for instance, would serve such a purpose. The focus of our study was to explore possible advantages of MECE and MECE + TV over SWE in the reconstruction of typical structures that have clinical relevance. Extending this study to test phantoms with other types of structures and elasticity values can provide additional insight into the potentials of the proposed method and more statistically meaningful results. It is also noted that the present study was based on and limited to shear elasticity reconstruction. Extension of MECE to viscoelasticity reconstruction can provide additional information about tissue mechanical properties. There are other aspects of MECE that need to be studied. Since MECE is a full-wave reconstruction, it is expected to be more computationally demanding than SWE. The reconstruction times for one data set with MECE are in the order of a few minutes compared to seconds in SWE. However, the algorithm performance can be significantly improved through simple steps such as pre-computation of element matrices, reduced-basis for dimensionality reduction, etc. We plan to do this in future versions of the method. Further work is needed to improve the computational efficiency of MECE.
In conclusion, based on this phantom study, the proposed MECE and MECE + TV elasticity reconstruction methods may offer better performance over the conventional SWE in terms of geometrical shape recovery, elasticity value estimation, and image contrast. With these advantages, the proposed method may have a significant impact on clinical application of elastography. The results of this study warrant further study of MECE and MECE + TV for in vivo applications.