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.

Curl‐based direct inversion

Direct inversion

Finite element method

Local frequency estimation

Multifrequency dual elastovisco

Magnetic resonance elastography

Traveling wave expansion.

INTRODUCTION

Tissue stiffness is often affected both directly and indirectly by disease, with common examples being fibrosis, tumor growth, and inflammation. As tissue becomes fibrotic, collagen density increases in the extracellular matrix, leading to a macroscale effect on stiffness.1, 2 During tumor growth, angiogenesis, increased cellular stiffness, and compaction of surrounding tissue can combine to alter material properties.3, 4 The significant changes in hardness often accompanying disease have historically led to palpation being used as an early tool for the detection of certain diseases. While still invaluable for initial screening, palpation is qualitative, limited to superficial regions of the body, and often most applicable in advanced stages of disease.

Elastography is an imaging technique that infers stiffness by examining a tissue's response to some excitation. This excitation, often externally applied, leads to deformation in the tissue, which can be measured by medical imaging methods such as magnetic resonance (MR) or ultrasound (US). These measured deformations can be used to reconstruct stiffness values by solving equations based on a chosen tissue model that relates these quantities. The calculated values of stiffness (typically the shear modulus) are similar to what is felt during palpation.5 However, elastography, besides being quantitative, has the added benefits of being able to examine deeper in the body, potentially detecting small changes in stiffness present in earlier stages of disease, and measuring additional viscous properties of the tissue.6 Furthermore, the shear modulus has been seen to vary more over various healthy and diseased tissue types than the bulk modulus or other imaging modalities such as T 1 relaxation.7

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 fibrosis24, 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 breast28 and liver tumors from surrounding healthy and fibrotic tissue29 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 sclerosis33, 34 (as opposed to US elastography, which has only been used to measure brain stiffness in surgical contexts35, 36.

Figure 1

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 al24 with permission from Elsevier)

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.

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

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 (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/m3, when modeling tissue. The shear modulus, G, is often written as G=G ′+i G ′ ′, where G ′ is the storage modulus and G ′ ′ is the loss modulus. Generally, Equation 2 can be inserted into Equation 1, giving although this is typically only done for compressible materials. Often, tissue is modeled as incompressible and Equation 2 becomes the divergence‐free condition on the displacements, In a typical forward problem, the material parameters are known and u and p are unknown. Boundary conditions are also necessary. For example, Dirichlet conditions such as may be specified, where are known displacements.

The weak (or variational) form of these equations is also noteworthy, as it is the weak form on which the finite‐element method (FEM) operates. The FEM is a popular numerical discretization choice for MRE because, as an integral approach, it reduces the order of differentiation, and thereby may help to reduce noise sensitivity. For example, the incompressible forward problem defined by Equations 1 and 4 has the weak form: find such that, for all , with the trial function space and the test function space . 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, , 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 , where i counts over the discrete elements of , typically one for each voxel. The goal of reconstruction methods is to solve for G, based on the measured , and, depending on the method, potentially include solutions for either λ or p.

The measured displacement fields, , 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 MRE39. 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 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 works38, 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 , where are calculated displacements, are measured displacements, 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 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 al46 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).

Figure 3

Comparisons of brain MRE data reconstructed using conjugate gradient, quasi‐Newton, and Gauss–Newton methods and compared with A, a T 2‐weighted image. B–D, G ′ reconstructions are shown using the compressible equations for the three methods in the order previously listed. E–G, G ′ reconstructions are also shown using the incompressible equations. H–M, The corresponding map of G ′ ′. Values are in units of kPa (©2017 IEEE. Reprinted, with permission, from Tan et al46

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 further48 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 described52 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 developed58 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). Doyley et al52 reported computation times of 1.5–4.5 hours by varying the subzone size (twelve 2.5‐GHz Xenon IV compute nodes running Linux with problem sizes greater than 104 FEM nodes). Perreard et al60 reported approximately 5 hour run times (eight AMD opteron processors running 2.5×105 element meshes with 100 global iterations.)

Table 2

Many in vivo and ex vivo human and animal studies using MRE have been published and a selection are listed in this table. The location in the table indicates the area of the study (row) and the type of reconstruction used (column). The first portion of the table includes studies on particular disease and treatment areas, while the second portion covers specific organs. The highlighted methods are as follows: LFE – local frequency estimation method; DI – direction inversion methods using high‐pass filtering or prescribed near‐incompressible values to address the compressional wave; CDI – direction inversion methods using the curl to address the compressional wave; MDEV – multifrequency dual elastovisco inversion method (Section 3.2.4); Subzone – iterative method as described in Section 3.1.2. The last category (Other) includes two studies using a global direct method with sparsity regularization120, 176 and one using TWE174

Area	LFE	DI	CDI	MDEV	Subzone	Other
Brain cancer			112	113
Breast cancer	114, 115	94	28, 116
Liver cancer	29		117
Prostate cancer	118		119			120
Liver fibrosis	24, 121, 122, 123, 124, 125		25
Alzheimer's			31		30
Parkinson's				32
Multiple sclerosis		33, 126	34	126	127
Ablation	128, 129
Liver response	130, 131
Tumor response		132	133
Brain	134	135, 136, 137, 138, 139, 140, 141	112, 142, 143, 144, 145	146, 147	148, 149, 150, 151
Breast	152
Heart	153		154
Kidney	155, 156, 157		158	159
Liver	160	136
Lung	161, 162, 163, 164	165
Muscle		166, 167	168, 169, 170, 171
Pancreas	172		173
Prostate	174	175				174, 176

Gauss–Newton methods

Many other reconstructions based on Gauss–Newton iterations have been published, in addition to the subzone technique. Kallel and Bertrand61 solved the static elasticity equations in 2D with Tikhonov regularization and near‐incompressible values of ν. Later, this approach was modified to utilize total variational regularization.62 Doyley et al63 and Miga64 followed similar approaches and used Levenburg–Marquardt. While the previous studies were oriented towards US, Khalil et al18 followed this approach for OCT with Tikhonov regularization.

Eskandari et al65 solved the time‐harmonic equations in a similar fashion using Levenburg–Marquardt. This implementation was later used in a study on bandpass sampling66 and was extended to employ mesh adaptation.67 Honarvar et al68 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 al69 introduced the traveling‐wave expansion (TWE) method for MRE. This method fits fundamental solutions of the Helmholtz equation to the measured displacement data and also assumes local homogeneity of stiffness. The solution chosen for fitting is a summation over a finite set of directions of complex exponential functions with a constant (and complex) wave number. The method works voxel by voxel and, for any wave number, a linear minimization is performed to calculate the unknown parameters in the solution function, thereby fitting the solution to some local subset of the data. The method is categorized as iterative here, because it performs an additional minimization over possible complex wave numbers, which 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 104 voxels in phantom studies,69 although no CPU information appears to be given.

Other

Samani et al70 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 al73 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 is the displacement data vector and 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, 75, 76 Eskandari et al76 suggested a quadratic programming technique, which constrains stiffness values to physically reasonable ranges. An adjoint weighted equation formulation was presented and demonstrated in multiple instances such as compressibility and anisotropy.77, 78 Zhang et al79 extended this approach 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 Maniatty40 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 al76 and the following heterogeneous methods.

Honarvar et al81 suggested a similar process to reconstruction to Park and Maniatty,40 but used sparsity regularization instead of Tikhonov for both pressure and stiffness. In Honarvar et al, the sparsity regularization is performed using the discrete cosine transform and masking in frequency space to limit the stiffness and pressure distributions to lower frequency contributions, effectively smoothing the result. As opposed to Tikhonov, sparsity reduces the size of the matrix system, but also globalizes the solution, filling the otherwise sparse and banded matrix. Removing these properties of the system should lead to a significant increase in computation time, unless only very few frequencies are included. On the other hand, filling the matrix could help to stabilize the solution.

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 Computation times of 15 minutes were reported for the curl form for 3D data and 35 minutes for the mixed form82 (computed on 15×15×15 overlapping subdomains of a 140×46×30 ROI with a MATLAB implementation, but CPU information was not reported). In Fovargue et al,38 computation times of approximately 50–90 seconds were reported for data sets with approximately 24 000–33 000 voxels for a mixed heterogeneous implementation with Tikhonov regularization (MATLAB with an Intel Xeon 4‐core, 8‐thread, 3.6 GHz processor). In Eskandari et al,76 compute times for the 2D heterogeneous direct method were compared with those of an iterative method65 and local method for varying image sizes (Figure 4).

Figure 4

Plot comparing the computational cost of various methods considering 2D data (measured on a quad‐core CPU with two 2.82‐GHz processors): a single iteration of an iterative method,65 a heterogeneous direct method, a local direct method, and the quadratic programming technique outlined in Eskandari et al.76 (©IOP Publishing. Reproduced with permission. All rights reserved. Eskandari et al. 2011)

Local frequency estimation

Manduca et al83 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 al85 investigated the use of Gauss filters in place of the lognormal quadrature filters used previously83 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 al88, 89 using a variational approach and Oliphant et al90, 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 al93 also presented spatio‐temporal directional filtering as a pre‐processing step and compared results from direct inversion and LFE algorithms. Sinkus et al94 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 al96 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 al41 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, 8‐thread, 3.6‐GHz processor).

Direct inversion methods have been used extensively in both research and clinical studies. McLaughlin et al97 used this approach to investigate filtering. Papazoglou et al98 investigated factors such as SNR and actuator orientation on the resulting elastogram and Riek et al99 investigated the frequency dependence of stiffness. 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 for the frequency dependence of stiffness and potentially longer scan times. Papazoglou et al100 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 al101 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 al102 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.

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 al101

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 amplitude103, 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 fits107 or simple peak‐to‐trough measurements.108, 109, 110, 111

Table 1

A summary of many of the reconstruction methods mentioned in this review. Here, they are presented in chronological order and categorized in terms of their approach to the inverse problem, whether they assume local stiffness homogeneity, and their consideration of the compressional component of the wave data. For the latter category, the following symbols are used: IG – ignore; ∇× – compute the curl of the wave data; NI – prescribe a near‐incompressible value; P – additionally solve for pressure in the mixed displacement‐pressure formulation; HPF – apply a high‐pass filter to the data; Gλ – solve for both the shear modulus, G, and the first Lamé parameter, λ

Method	Inverse problem approach	Local stiffness homogeneity	Compressional component
Kallel et al61	Iterative	No	NI
Manduca et al83	Direct	Yes	IG
Romano et al88	Direct	Yes	IG
Oliphant et al90	Direct	Yes	IG
Van Houten et al47	Iterative	No	NI
Doyley et al63	Iterative	No	NI
Fu et al73	Iterative	Yes	NI
Sinkus et al94	Direct	Yes	NI
Manduca et al92	Direct	Yes	HPF
Samani et al70	Iterative	No	NI
Van Houten et al49	Iterative	No	NI
Miga et al64	Iterative	No	NI
Oberai et al44	Iterative	No	P
Zhu et al74	Direct	No	NI
Doyley et al52	Iterative	No	Gλ
Sinkus et al95	Direct	Yes	∇×
Park and Maniatty40	Direct	No	P
Eskandari et al65	Iterative	No	NI
Guo et al75	Direct	No	NI
Baghani et al69	Iterative	Yes	NI
Eskandari et al76	Direct	No	NI
Okamoto et al96	Direct	Yes	IG
Honarvar et al81	Direct	No	P
Papazoglou et al100	Direct	Yes	HPF
Zhang et al79	Direct	No	P
Honarvar et al82	Direct	No	∇×
McGarry et al51	Iterative	No	NI
Hirsch et al101	Direct	Yes	∇×
Connesson et al41	Direct	Yes	Gλ
Honarvar et al68	Iterative	No	P
Tzschätzsch et al102	Direct	Yes	HPF
Tan et al46	Iterative	No	P
Fovargue et al38	Direct	Yes	P

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 other researchers, such as the LFE method in the MRE Wave software (Mayo Clinic, Rochester, MN, USA) and Resoundant software (Rochester, MN, USA), as well as the MDEV, k‐MDEV, and LFE methods on a publicly accessible website.177 Several possible reasons for this disparity include low computation time, ease of use, and guaranteed solutions.

Local direct methods are by far the least computationally expensive; however, few publications on these methods38, 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 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 MRE46, 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.