Stiffness reconstruction methods for MR elastography

Assessment of tissue stiffness is desirable for clinicians and researchers, as it is well established that pathophysiological mechanisms often alter the structural properties of tissue. Magnetic resonance elastography (MRE) provides an avenue for measuring tissue stiffness and has a long history of clinical application, including staging liver fibrosis and stratifying breast cancer malignancy. A vital component of MRE consists of the reconstruction algorithms used to derive stiffness from wave‐motion images by solving inverse problems. A large range of reconstruction methods have been presented in the literature, with differing computational expense, required user input, underlying physical assumptions, and techniques for numerical evaluation. These differences, in turn, have led to varying accuracy, robustness, and ease of use. While most reconstruction techniques have been validated against in silico or in vitro phantoms, performance with real data is often more challenging, stressing the robustness and assumptions of these algorithms. This article reviews many current MRE reconstruction methods and discusses the aforementioned differences. The material assumptions underlying the methods are developed and various approaches for noise reduction, regularization, and numerical discretization are discussed. Reconstruction methods are categorized by inversion type, underlying assumptions, and their use in human and animal studies. Future directions, such as alternative material assumptions, are also discussed.

Anatomical, wave, and stiffness images of normal and fibrotic liver. While little difference is apparent in the anatomical images, large differences are seen in the wave data and reconstructed stiffness. Also clearly demonstrated is the relationship between increased wavelength and increased stiffness (Reprinted from Yin et al 24 with permission from Elsevier) FIGURE 2 An illustration of a typical process of harmonic MRE. A, The patient or object is subjected to a single frequency vibration while in the MR scanner. Here, data collected from MRE performed on a CIRS 049 elastography phantom (CIRS Inc., Norfolk, VA, USA) are shown. B, During phase contrast imaging, the complex-valued k-space data are collected and constructed into a magnitude image and 3D displacement field images at several time points along the cycle of the vibration. C, The real-valued time-dependent displacement data are converted by Fourier transform to complex-valued steady-state wave data. D, A reconstruction algorithm solves an inverse problem to compute the complex shear modulus. 38 The shape and differing stiffness values of the inclusions in the phantom are readily apparent in the reconstruction, whereas in the MR magnitude image in B the inclusions appear very similar. The imaginary part of the shear modulus in this phantom should be zero throughout. This appears to be the case with this reconstructed image, besides the small spike in the stiffest inclusion. This stiff region corresponds to a long wavelength and this highlights an important balance in MRE. For longer wavelengths, noise becomes a more dominant factor during calculations, for some reconstruction methods, potentially leading to more error in the resulting stiffness The first experiments in elastography began with US imaging. [8][9][10][11] The low cost and wide availability of US continues to make this a popular elastography imaging method. 12 Later, MR methods for elasticity imaging were developed, [13][14][15] with advantages including a full 3D image stack of the tissue displacements and large, movable, and deeper fields of view. 7,12 For a history of US and MR elastography (MRE), see Sarvazyan et al. 16 Elastography has also been performed with optical coherence tomography (OCT), where advantages include high resolution and smaller scale imaging, [17][18][19][20][21] X-ray, 22 and mechanical imaging, 23 which measures surface stress patterns to infer internal structure.
MRE has demonstrated success in many areas. Investigations into liver have revealed that MRE is capable of staging liver fibrosis 24,25 (Figure 1) with greater separation than other modalities, including US elastography and MR diffusion-weighted imaging. 26,27 MRE has also seen success in differentiating benign and malignant tumors in breast 28 and liver tumors from surrounding healthy and fibrotic tissue 29 by their viscoelastic characteristics. MRE is able to image waves in the brain, leading to research on stiffness changes due to Alzheimer's, 30,31 Parkinson's, 32 and multiple sclerosis 33,34 (as opposed to US elastography, which has only been used to measure brain stiffness in surgical contexts 35,36 ).
Within MRE there are many differing implementations, where considerations include the type of tissue excitation, the hardware required to produce the excitation, acquisition methods, and reconstruction methods. Three common tissue excitation techniques across elastography are quasi-static deformations, transient waves due to singular pulses, and harmonic wave motion due to a single frequency vibration; however, the most common technique in MRE is harmonic. 37 The process of harmonic MRE ( Figure 2) works by inducing periodically steady waves in a region of interest (ROI) in a patient or object, performing phase-contrast MR imaging, and, from the complex-valued k-space data, constructing displacement data representing wave motion within the ROI. Finally, a reconstruction method is applied to calculate the underlying stiffness that modulated the measured wave behavior. Therefore, elastography reconstruction is considered an inverse problem. Typically, due to the small displacements involved, the linear elasticity or linear viscoelasticity equations are used to relate the displacements and stiffness, and hence it is these equations that are inverted. The resulting image of elasticity is usually of the same resolution as the displacement images and is sometimes called an elastogram.
MRE reconstruction algorithms have been implemented in many ways, with different underlying assumptions and employing varying processes for the inverse solution. Understanding these differences is important, as they effect uniqueness, accuracy, robustness, ease of use, and required computation time. This article discusses these effects while comparing and summarizing MRE reconstruction methods and associated studies currently present in the literature. The theory underlying MRE reconstruction is presented and the applicability of MRE in various human and animal studies is described.
In Section 2, theory and background that is common across many MRE reconstruction methods is presented. Section 3 contains descriptions and categorizations of various methods and includes some mathematical underpinnings and references to reconstruction usage in other studies.
A discussion and conclusion are provided in Section 4, highlighting the scientific and clinical value of MRE as well as some of the key challenges in reconstruction.

THEORY
In harmonic MRE, a transducer is placed against a patient and produces single-frequency waves (typically 30-90 Hz) that move through an ROI. Due to the small displacements associated with these waves, the linear elasticity or viscoelasticity equations are typically assumed to govern the tissue deformation. An isotropic assumption is also typically made of tissue during MRE, which is usually considered reasonable depending on the tissue type. The isotropic linear elastic equations contain three material parameters. The first is the density, , and common pairs for the remaining two include the following: (i) the two Lamé parameters, and G (the latter is also called the shear modulus); and (ii) the Young's modulus, E, and Poisson's ratio, . As materials approach incompressibility, the following occur: → 0.5, → ∞, and E → 3G. For tissue that is nearly incompressible, knowing the shear modulus is therefore essentially equivalent to knowing the Young's modulus.
The wave displacements are constructed from complex-valued k-space MRI data. These data are captured at several time points over the period of the wave (four or eight time points is common). Occasionally, even at this stage, qualitative conclusions can be made regarding the underlying stiffness, as changes in wavelength are sometimes visible and these correspond to changes in stiffness (in the absence of boundary effects). For quantitative results, however, reconstruction is used. Typically, due to the time harmonic excitation, dynamic vibrational analysis is used, which models the displacements and stiffness parameters as complex-valued entities. The equations, in this form, then also capture effects of viscoelasticity.
A discrete Fourier transform converts the time-dependent real displacements, u r (x, t), into a steady-state complex-valued field, u(x), related by , and reconstruction methods calculate complex values for stiffness, G.
The isotropic linear viscoelasticity equations are where u, p, G and are complex valued. Additionally, p is pressure and is angular frequency. The density, , is usually assumed to be that of water, ≈ 1000 kg/m 3 , when modeling tissue. The shear modulus, G, is often written as G = G ′ + iG ′′ , where G ′ is the storage modulus and G ′′ is the loss modulus. Generally, Equation 2 can be inserted into Equation 1, giving Here, w and q are test functions, Ω is the spatial domain, Γ 1 denotes the portion of the boundary of Ω where known displacements,ū, are imposed, and Γ 2 denotes the portion where known tractions,T, are used.
However, in MRE an inverse problem is solved, where u, , and are known and G, , and p are unknown. Boundary conditions may also need to be reconsidered depending on the inverse formulation. Furthermore, only an estimate of u is measured from MR; this is a discrete quantity that contains error (e.g. noise or artifacts), herein denotedū i , where i counts over the discrete elements ofū i , typically one for each voxel. The goal of reconstruction methods is to solve for G, based on the measuredū i , and, depending on the method, potentially include solutions for either or p.
The measured displacement fields,ū i , inherently contain noise, as MRI produces images with finite signal-to-noise ratio (SNR). A common and straightforward method to reduce noise is to apply a local Gaussian filter directly to the displacement data. Another common approach is local polynomial fitting, also called Savitzky-Golay filtering. From this method, smooth derivatives of the displacement fields can also be computed by using the derivatives of the polynomials calculated during the filtering. These local spatially based filters may have issues near discontinuities and boundaries, but these errors will be localized. Frequency-based low-pass filtering is another option, but this may have issues in oddly shaped domains, and, being a global operation, errors may propagate throughout the domain. More advanced denoising techniques have also been applied in MRE 39 .
Noise, in part, limits the resolution of MRE for many reconstruction methods, as high-resolution images would lead to many pixels per wavelength and noise would then dominate many calculations, instead of the curvature of the waves themselves. The other part of this balance is a limit to decreasing the wavelength, as this corresponds to increasing the frequency, which will increase attenuation.
Both compressional and shear waves are present in tissue during MRE and it is the combined displacements from both that are collected by imaging. The assumption of linear viscoelasticity accounts for both types of wave; however, it is the shear modulus, and therefore the shear waves, that elastography is oriented towards. In part, this is due to the shear-wave speed varying more over different tissue types than the compressional wave speed, but also to the fact that the compressional waves are much faster, which makes balancing considerations like applied frequency and imaging constraints more difficult. 6 However, since the compressional wave is represented in the imaged displacements, it must somehow be considered during reconstruction. In the equations, the compressional component corresponds to Equation 2 and the ∇p term in Equation 1 or ∇ (λ∇ · u) in Equation 3.
Amongst reconstruction methods, there are five common approaches for considering the compressional wave: (1) ignore this term, i.e. assume the gradient of pressure is negligible; (2) apply a discretized curl operator to the displacement data and reconstruct using this curl field, assuming it to be divergence-free; (3) prescribe one of the parameters, most often by recasting Equation 3 in terms of and E, assuming a near-incompressible value for , and solving for E; (4) apply a high-pass filter aimed at removing the long-wavelength compressional wave; and (5) assume incompressibility and solve Equation 1 for p in addition to G.
There are strengths and drawbacks to each of these approaches. Neglecting the pressure term has been shown to lead to errors in the reconstructed stiffness, 40 and, while this approach is mentioned as a possibility in several early MRE manuscripts, it is less common in contemporary methods. Applying the curl and attempting to limit hydrostatic effects increases the order of derivatives applied to the displacement field. Since the displacements contain noise and differentiation amplifies noise, this presumably increases the sensitivity of the reconstruction. Using an accurate near-incompressible value of typically leads to numerical instability, so usually values in the range 0.49-0.499 are chosen. This causes the errors in the reconstructed to be of several orders of magnitude. Reconstructed G may be accurate, however there may be errors in regions of strong mode conversion.
When applying a high-pass filter, it is not always the case that the compressional and shear components are well separated in frequency space.
The shape of the domain and areas of mode conversion can lead to a lack of separation. Finally, solving for p increases the number of unknowns, thereby demanding more from the available data. The pressure in some cases may also require some regularization. Two works 38,41 have suggested using specialized virtual fields and FEM test functions, respectively, to remove the compressional component from the equations. This seems to be most similar to including the second unknown (bulk modulus or pressure) in the solution, but could be considered as its own category.
While the preceding discussion has focused entirely on theory assuming harmonic excitation, as this is most common in MRE, some of the methods described in Section 3 assume quasi-static deformation. Both of these approaches assume a steady state and so techniques from one category may be applicable to the other.

METHODS
Reconstruction methods for MRE may be categorized in a variety of ways, but probably the largest division regards the two approaches to solving the inverse problem: iterative and direct. In other literature, iterative methods are sometimes called nonlinear inversion (NLI) methods.
Iterative methods solve a nonlinear constrained minimization problem, that is they attempt to minimize the difference between the measured wave field and simulated wave fields found from solving forward problems. During each iteration, a forward problem is solved using a set of stiffness variables calculated from a previous iteration using minimization update techniques. Upon reaching a minimum, or satisfying some convergence criteria, the current stiffness variables are considered to be the solution to the inverse problem. Iterative methods, therefore, are strongly dependent on model assumptions such as initial stiffness values and boundary conditions, as they only consider displacement fields that satisfy the governing equations. Features in the landscape of the stiffness variable space, such as multiple minima or shallow minima, may lead to a lack of uniqueness, slow convergence, or divergence. Regularization techniques are often employed to mitigate these issues, as well as decreasing the likelihood of inaccurate or non-physical results. These come with the drawback of additional parameters, on which the results will also depend. Iterative methods are currently computationally expensive, due to the many matrix solutions required for computing both the forward solutions and the stiffness variable updates; however, many potential advances may address this. Despite these concerns, iterative methods are a useful approach, due to their relative independence from data quality and their potential for increased accuracy.
In contrast, direct methods, while having many benefits, are more dependent on data quality. They solve a linear minimization problem, by assuming that the displacements from the measured wave field are sufficiently accurate to be inserted into the governing equations, leaving the linearly dependent stiffness variables (and potentially pressure) as the unknowns available to minimize the system. After computing any required displacement derivatives by standard numerical techniques, the system is solved directly for the unknowns. Typically, the resulting matrix system is overdetermined and so is solved by least-squares matrix inversion. Due to the single solution, these methods require much less computational effort than iterative methods. Other advantages include fewer method parameters and the guaranteed existence and uniqueness of a solution. The solution itself is, however, more sensitive to data quality and any processing parameters, compared with iterative methods.

Iterative methods
The goal of minimization in iterative methods is to find subject to (û, p, G) = 0, whereû are calculated displacements,ū are measured displacements, ‖·‖ X denotes some norm, and the operator, , is presumably based on the equations of motion presented in Section 2.
The measured displacement data are often used in two ways within iterative methods. The first is to define the optimal configuration of the displacements calculated from the forward problem, as in Equation 7. The second is to specify the boundary conditions for the forward problem. Along with boundary conditions, an initial distribution of stiffness must also be chosen, at which point a forward problem may be solved for displacements via numerical discretization techniques, such as FEM. From here, some optimization algorithm is followed in order to update the stiffness.
This typically requires computing changes in the solution due to changes in stiffness, which, in turn, requires solving a similar forward problem for each discrete stiffness parameter. The computed stiffness update is applied and a new forward problem is solved with this updated distribution.
This process repeats until some stopping criterion is reached. Common minimization approaches include gradient descent, Gauss-Newton, and Newton's method. In this order, these methods are increasing in the amount of computation required per iteration but decreasing in the number of iterations required in theory. It seems that the most common in MRE reconstruction is Gauss-Newton, as this does not require construction of the second-order derivatives for Newton's method but offers a more accurate update than gradient descent.
One straightforward and common way to solve Equation 7 is least-squares, which attempts to minimize the function where G is a vector containing the material parameters, i counts over voxel locations, j counts over the three components of the displacements, and * denotes the complex conjugate. Derivatives of F with respect to each material parameter are computed and an update formula is constructed, typically of the form where ΔG n is found by solving (for a Gauss-Newton approach) Here, r is the residual vector and J contains the derivatives of F.
A common modification to the Gauss-Newton optimization algorithm is Levenberg-Marquardt. 42 This approach weights the approximate Hessian, J T J, towards a gradient descent system by an adaptively updated weighting parameter. This may be thought of as regularization of the approximate Hessian, but, importantly, this strategy does not change the stiffness variable landscape, only the approach towards a minimum. Considering the form presented in Equation 10, the gradient-descent method simply replaces J on the left-hand side with the identity and typically a multiplier is used to limit the step size. A full Newton's method approach, which would replace J T J with the Hessian matrix, has not yet been applied to elastography. Other options include quasi-Newton methods, such as BFGS, which approximate the Hessian as the iterations proceed.
Regularization of the stiffness parameters is also commonly employed, often by adding some measure of the derivatives of the stiffness to the minimization function, thereby encouraging a smoother result. Tikhonov regularization, for example, involves adding the L 2 norm of the Laplacian of the stiffness. 43 These techniques usually require this additional regularization term to include a weighting parameter to weight this portion of the minimization against the norm of the displacement differences. This type of regularization changes the stiffness variable landscape, favoring smoother solutions and ideally improving the likelihood of a clear minimum.

Gradient-descent methods
The application of gradient-descent methods to elastography has been investigated by Oberai et al, 44,45 in which an adjoint operator is used to modify the linear elasticity equations to realize the descent calculation. This is more efficient than standard gradient-descent methods, as it requires only two solutions for each iteration. The mixed displacement-pressure formulation of the static deformation equations are solved via FEM in 2D and the minimization is regularized via a Tikhonov approach. Tan et al 46 introduced a generalized adjoint method applicable to poroelasticity and also provided comparisons between conjugate-gradient, quasi-Newton, and Gauss-Newton methods by computational complexity and by results across various data sets, including phantom and brain ( Figure 3).

Subzone
A particularly prominent iterative reconstruction method, and one usually oriented towards harmonic MRE, is the subzone technique, proposed by Van Houten et al. 47 In later publications, this method was evaluated further 48 and extended to 3D. 49,50 This method has typically been implemented with Gauss-Newton, but also with quasi-Newton methods, 51 and usually employs Levenberg-Marquardt. This approach uses a FEM discretization of the equations and performs a minimization procedure in many small overlapping domains, or subzones, that span the global domain. Both local and global convergence are usually considered. While specifics depend on the particular implementation, typically some minimization iterations are performed in the subzones and, at certain intervals or after local convergence, the subzone domains are changed and values are exchanged between them. Another round of minimization within subzones is performed and the process repeats until a global convergence criterion is satisfied.
In various implementations of this method, either is set to a near-incompressible value and E is solved for, both and E are solved for, or both and G are solved for.
Many updates and extensions to this method have been published. A parallelized 3D implementation was described 52 and parameter choices for this implementation were analyzed later. 53 Density was added as an unknown parameter, 54 a Rayleigh-damped model has been considered, 55 and poroelasticity has been investigated. 56,57 A multiresolution approach was developed 58 and soft prior regularization was presented, which predefines areas of homogeneous material parameters. 51 The method has been used in many studies, including those on SNR, 59 frequency and direction dependence in phantoms, 60 and various other in vivo and ex vivo data ( Table 2) ing the subzone size (twelve 2.5-GHz Xenon IV compute nodes running Linux with problem sizes greater than 10 4 FEM nodes). Perreard et al 60 reported approximately 5 hour run times (eight AMD opteron processors running 2.5 × 10 5 element meshes with 100 global iterations.)

Gauss-Newton methods
Many other reconstructions based on Gauss-Newton iterations have been published, in addition to the subzone technique. Kallel 66 and was extended to employ mesh adaptation. 67 Honarvar et al 68 solved the mixed displacement-pressure and time-harmonic form using Gauss-Newton and Levenburg-Marquardt. This implementation was also used for comparison with heterogeneous direct methods. 68

Traveling-wave expansion
Baghani et al 69  is achieved by a Newton-type iteration, so, at each voxel, the method computes the wave number that minimizes the error between the data and a local fit of the data. The wave number for the compressional component is prescribed. Although an iterative solution is required at every voxel, the iterations only involve two parameters and, when fitting, the system sizes are small due to the local constructions. This keeps the computation time low and times of less than one minute are reported for image sizes of the order of 10 4 voxels in phantom studies, 69 although no CPU information appears to be given.

Other
Samani et al 70 introduced a method that uses the linear elasticity equations themselves as an update formula for iterating towards the solution.
This method uses FEM and considers the quasi-static 3D equations with set to a near-incompressible value. A study published later used and extended this method. 71 A genetic algorithm approach was presented by Zhang et al, 72 which also uses FEM and a near-incompressible . Fu et al 73 considered local stiffness homogeneity within an iterative approach. This method performs a minimization procedure at each voxel by considering increased and decreased stiffness and finding the midpoint. This method assumes 2D and a near-incompressible . As with the TWE approach, the homogeneity assumption significantly reduces computational time.

Direct methods
Direct methods take advantage of the abundance of data acquired in MRE to solve directly for the stiffness. This leads to a simpler approach than iterative methods, usually requiring fewer parameters and much less computational effort. If stiffness is allowed to vary spatially, then regularization is sometimes still required to stabilize the solution. These heterogeneous direct methods typically solve over the entirety (or large portions) of the domain and attempt to capture the true heterogeneous nature of the stiffness. In contrast, local direct methods assume local homogeneity of stiffness and perform independent solutions on many small subdomains that combine to span the full domain.
While the accuracy of local methods will suffer in regions with significant stiffness heterogeneity, the homogeneity assumption simplifies and stabilizes the methods as well as decreasing computational expense. In the simplest case of local direct inversion, and where the compressional component is assumed to have been removed, the stiffness at voxel i could be found by least squares: whereū i is the displacement data vector and (Δu) i is some approximation of the Laplacian of the data at voxel i.

Heterogeneous direct methods
Heterogeneous direct methods consider stiffness to vary spatially and invert the governing equations directly with the proper derivatives affecting the stiffness. These methods are less computationally expensive than iterative methods, as only a single solution is required, but are more expensive than local direct methods, as the matrix system is much larger. Allowing for stiffness heterogeneity may lead to more instability and sensitivity in the method, and therefore these methods sometimes require regularization. As with iterative methods, a common regularization approach is Tikhonov, where, again, a measure of the Laplacian of the stiffness is added to the minimization. If pressure is additionally solved for, then this may also be regularized.
Several methods consider 2D quasi-static elasticity with near-incompressible values of and discretization via FEM. 74 to the time-harmonic case in 2D with a mixed formulation and total variation regularization. Effects of noise on this method were presented later. 80 Park and Maniatty 40 solved the mixed displacement-pressure form of the harmonic linear elasticity equations. In this approach, FEM is used with Tikhonov-type regularization applied to the pressure, whereas the stiffness is not regularized. The measured displacement data are filtered by solving for a new displacement field that minimizes three terms: the difference compared with the measured displacements, the divergence of the field, and the Laplacian of the field. Between the pressure regularization and displacement adjustment, at least five regularization parameters are introduced. In the FEM solution, all equations associated with boundary nodes are removed, which removes the unknown boundary terms in Equation 5. This technique was also used in Eskandari et al 76  A similar method was later introduced in which the curl operator is applied to the equations to remove the pressure term. 82 In the finite-element process, integration by parts is applied for an additional time to what is standard in order to remove the extra derivative introduced by the curl.
Then, in order to remove the boundary conditions, special basis functions are utilized. Similar results were reported between the mixed and curl forms, although, since removing the pressure term reduces the number of unknowns, the curl form is computationally faster. 82 The mixed form was later compared with an iterative method. 68  In Eskandari et al, 76 compute times for the 2D heterogeneous direct method were compared with those of an iterative method 65 and local method for varying image sizes (Figure 4).

Local frequency estimation
Manduca et al 83 presented a local direct method called the local frequency estimation (LFE) method, which is based on particular frequency banks and filtering described in Knutsson et al. 84 The stiffness is then calculated from the frequency estimate by assuming homogeneity and no attenuation (meaning only the real part of the complex stiffness modulus is provided). Although the compressional component is not accounted for directly, its effect may be minimized by limiting the frequency range during filtering. Braun et al 85 investigated the use of Gauss filters in place of the lognormal quadrature filters used previously 83 and principal frequency estimation was proposed by McGee et al. 86 LFE has also been used to investigate the effect of off-frequency sampling, 87 as well as in many ex vivo and in vivo animal and human studies ( Table 2).

Direct inversion
Direct inversion methods (sometimes denoted DI) perform a pixel-wise inversion of the homogeneous time-harmonic linear viscoelasticity equations. This was first suggested in Romano et al 88,89 using a variational approach and Oliphant et al 90,91 using polynomial fitting to compute the necessary derivatives. In these first implementations, the compressional component is sometimes ignored. Manduca et al, proposed high-pass filtering of the wave data in order to limit the low-frequency compressional component. 92,93 Manduca et al 93 also presented spatio-temporal directional filtering as a pre-processing step and compared results from direct inversion and LFE algorithms. Sinkus et al 94 followed a similar procedure to find tensor components of stiffness with a near-incompressible value for . Another common approach among direct inversion implementations is to apply a discrete curl operator to the displacement data and reconstruct using the curl field. 92,95 Okamoto et al 96 proposed a finite-difference-based direct inversion using total least squares, which may provide more accuracy when there is uncertainty in the independent variable. This work also compared MRE with dynamic shear testing in soft gels. Connesson et al 41 implemented a virtual fields method approach, wherein the virtual field is constructed to remove the bulk term and reduce dependence on noise. Fovargue et al. 38 proposed a FEM-based inversion that uses a specific basis for the FEM test functions to remove the pressure term. Computation times for local direct methods are also supplied in Fovargue et al, 38 where times ranging from 3-9 seconds were reported (MATLAB with an Intel Xeon 4-core,  Table 2 contains many examples of DI methods applied to human and animal studies, where, to provide some division, curl-based approaches are listed separately from non-curl-based approaches.

Multifrequency direct inversion
Combining multiple-wave data sets corresponding to the application of different vibrational frequencies is one potential avenue for increasing the stability and reliability of MRE, by essentially increasing the amount of data per unknown. This comes with the disadvantages of needing to account FIGURE 5 Diagram showing the process used in the MDEV inversion. The curl components (columns A and E) are used to find the magnitude and phase of G for each multifrequency data set (columns B and D), which is averaged in column C by comparison with a T 2 -weighted image (©2013 Wiley Periodicals, Inc. Reproduced, with permission, from Hirsch et al 101 ) for the frequency dependence of stiffness and potentially longer scan times. Papazoglou et al 100 proposed a local direct inversion method using high-pass filtering to limit the compressional component and utilizing multiple data sets representing different frequencies. The stiffness is computed in accordance with a springpot power law after first computing the exponent parameter. Hirsch et al 101 proposed a similar approach called multifrequency dual elastovisco (MDEV) inversion. The phase and magnitude of the complex shear modulus are computed separately and using averages over the multifrequency data sets. This method also computes and uses the curl of the wave data. A diagram illustrating the process is shown in Figure 5. MDEV is also used as the inversion algorithm in the elastography pipeline presented by Barnhill et al,39 as well as some studies on human data ( Table 2). Tzschatzsch et al 102 introduced a 2D multifrequency approach called k-MDEV, which uses filters in frequency space to construct a set of plane-wave approximations in multiple directions (these filters also work to limit the compressional component). Estimates of the stiffness and attenuation are found from plane-wave assumptions and by averaging over all frequencies and plane-wave directions. Multifrequency inversion is also applicable to other reconstruction methods, not only local direct inversion, and this is an area of active research.

DISCUSSION AND CONCLUSION
The methods presented in Section 3 are also summarized and categorized in Table 1, based on approaches to the inverse solution, local homogeneity assumption, and compressional component of the wave data. Some of the methods have been successfully applied to ex vivo and in vivo human and animal data. A summary of a selection of these studies is provided in Table 2. MRE has covered many organs, each with their own challenges as well as specific disease and treatment response areas. The liver fibrosis and brain sections of Table 2, for example, contain many studies, as these have shown much promise with MRE. More challenging areas include heart and skeletal muscle. Heart data often do not contain enough pixels to apply traditional reconstruction techniques, so, in addition to the one study listed in Table 2, other studies have measured relative changes in stiffness by examining changes in shear-wave amplitude [103][104][105] and by utilizing a spherical shell model. 106 Due to the strong anisotropy present in muscle, some methods, after carefully considering the wave propagation direction, apply damped sine-wave fits 107 or simple peak-to-trough measurements. [108][109][110][111] Although the iterative subzone, TWE, and global direct with sparsity regularization methods have been used in a few studies, local direct methods seem to dominate Table 2. Also, these methods, particularly LFE, have been implemented more often in software made available to clinicians and Local direct methods are by far the least computationally expensive; however, few publications on these methods 38,76 provide computation time estimates. This is presumably because they are too short to be of concern. Even for larger 3D data sets, the reconstruction time is measured in seconds (this is not including steps such as phase unwrapping; only stiffness reconstruction from properly processed displacement data). In contrast, computation times of heterogeneous direct methods are typically measured in minutes. These methods sometimes rely on regularization parameters being optimized, so may need to be run multiple times. This may also be true of the even more expensive iterative methods. The computation times of these methods, however, may vary considerably, as they depend strongly on implementation, domain size, and parameters. Times measured in hours have been reported, but with increasing computing power, domain-size reduction, and algorithmic advances, those times may be significantly reduced. Overall, few direct comparisons of computational time between methods have been reported and the comparisons made here should be taken as mostly anecdotal, as the reported times span years and are performed on different implementations with different CPU specifications.
Furthermore, local direct methods are usually not dependent on parameters, besides those controlling the initial smoothing of the wave field.
While some options may be available to researchers, these methods are typically able to be used out of the box. Methods that consider stiffness heterogeneity, whether direct or iterative, usually require regularization terms to restrain the stiffness. These terms are weighted against the original equations, but the weights of these parameters are difficult to know a priori and can affect the resulting elastogram. Iterative methods, being  112 113 Breast cancer 114,115 94 28,116 Liver cancer 29 117 Prostate cancer 118 119 120 Liver fibrosis 24,121-125 25 Alzheimer's 31 30 Parkinson's 32 Multiple sclerosis 33 more complicated, are additionally dependent on parameters other than stiffness regularization, such as subdomain size. While the parameter governing the Levenberg-Marquardt procedure is updated during the iterations, the initial choice may still have an effect. Results and comparisons from iterative and heterogeneous methods can therefore be harder to justify if significant parameter adjustment has occurred.
Reconstruction methods are sensitive to inconsistencies in the wave data, such as noise, but also imaging artifacts or other wave behavior that deviates from material assumptions. It is the regularization mentioned above that decreases sensitivity to these phenomena and allows methods to attain a level of well-posedness. Local direct methods again excel here. This is due to both the homogeneity assumption, which increases stability and robustness, and the direct approach, which guarantees a unique solution.
Although local direct methods dominate the studies in Table 2, there has been much research devoted to other approaches, as seen in Table 1.
This is presumably because of the potential for increased accuracy from heterogeneous and iterative methods. Heterogeneous approaches correctly capture changes in stiffness not only at large discontinuities but also in continuously varying regions. Local direct methods are inaccurate in these regions, to varying degrees, and changes in one component of the complex shear modulus can lead to inaccuracies in both. Furthermore, local methods are not guaranteed to predict an average of the underlying varying stiffness, but instead may calculate incorrect values, which is particularly evident near discontinuities. Iterative methods, additionally, are less sensitive to noise in the data, at least in theory. This also leads to increased accuracy, as noise tends to decrease computed stiffness in direct methods, because noise is seen as high-frequency waves, which would signal a lower underlying stiffness. Varying data quality is common too, as attenuation of the waves leads to decreasing SNR as distance from the transducer increases.
If iterative and heterogeneous methods can increase ease of use and robustness, they may begin to compete with local direct methods in terms of clinical viability. This may occur through approaches for quality assessment and confidence estimates or through some automatic adjustment of parameters. Perhaps too often, new reconstruction methods, or updates to existing ones, are validated only on in silico and phantom wave data, instead of being tested on more difficult anatomical data. Of course, these broader comparisons require well-established stiffness results or publicly available MRE data sets. Another option is the development of in silico data sets that properly capture the added difficulty inherent in anatomical data.
MRE reconstruction may also be extended by considering other material assumptions. Anisotropy is a common feature of tissue and so has been included in some reconstructions. 170,171,[178][179][180][181] Associated studies include investigations into the validity of anisotropic material assumptions, 182 results of applying isotropic reconstructions to anisotropic data, 60,150 experimental requirements, 183 and hardware to produce adequate wave data. 184 Poroelastic material assumptions have also been studied in MRE 46,56,57 as well as US elastography. 185 Nonlinear, usually neo-Hookean, constitutive laws have been investigated in the context of static elastography. [186][187][188][189] Many research directions within MRE are progressing as it continues to strengthen as a diagnostic method. It is becoming more widely used in tangential research and in the clinic through MRE packages available for MR systems. Success has been shown in staging liver fibrosis and potentially in stratification of tumor malignancy and various brain diseases. As precision and accuracy improve, it has the additional potential to replace invasive biopsy procedures for diagnoses. Stiffness reconstruction algorithms have played a vital role in positioning MRE at the forefront of medical imaging research through advances in robustness, accuracy, ease of use, and computational cost. Further improvements in these areas can allow MRE usage to become more widespread by increasing reliability and confidence in the approach.