Recovering vector displacement estimates in quasistatic elastography using sparse relaxation of the momentum equation

Abstract

We consider the problem of estimating the 2D vector displacement field in a heterogeneous elastic solid deforming under plane stress conditions. The problem is motivated by applications in quasistatic elastography. From precise and accurate measurements of one component of the 2D vector displacement field and very limited information of the second component, the method reconstructs the second component quite accurately. No a priori knowledge of the heterogeneous distribution of material properties is required. This method relies on using a special form of the momentum equations to filter ultrasound displacement measurements to produce more precise estimates. We verify the method with applications to simulated displacement data. We validate the method with applications to displacement data measured from a tissue mimicking phantom, and in-vivo data; significant improvements are noticed in the filtered displacements recovered from all the tests. In verification studies, error in lateral displacement estimates decreased from about 50% to about 2%, and strain error decreased from more than 250% to below 2%.

Keywords:

• Elastography
• biomechanical imaging
• strain estimation
• lateral displacement
• displacement recovery
• elasticity
• full-field

AMS Subject Classifications:

• 35R30
• 74G75
• 74B05
• 74L15
• 35Q74

1. Introduction

Ultrasound elasticity imaging (UEI), or elastography, is a rapidly growing field, primarily due to the increasing interest in the non-invasive quantification of the mechanical properties of soft tissues.[1–6] A crucial step in UEI is the estimation of a tissue’s motion. This motion is typically generated by a quasistatic compression, and the deformation of the tissue is measured with ultrasound. The measured deformation and an appropriate mathematical model can then be used to infer the mechanical properties of the soft tissue.

Ultrasound measurements of tissue deformation are typically much more precise in the ‘axial’ direction (i.e. the direction of the ultrasound beam), than in the so-called ‘lateral’ direction, the image direction transverse to the direction of the ultrasound beam.

The origin of the difference in precision is due to the anisotropy of the ultrasound point spread function (PSF). The width of the peak of the radio frequency (RF) PSF in the axial direction (the direction of sound propagation) is roughly an order of magnitude smaller than that in the lateral direction.[7]

With only one component of the displacement field measured precisely, it is impossible to compute precise estimates of the full strain tensor field, which would otherwise be useful information to researchers who do strain imaging. Furthermore, one goal of measuring the displacement field is as input to an inverse problem to estimate material property distributions.[6,8] Computing the material properties from only one component of the displacement field is computationally intensive.[6,8] In that situation, an iterative inversion approach is required, which may take hundreds to thousands of iterations to converge. Such an iterative inversion code usually needs displacement boundary conditions to drive the forward problem, and using the noisy measurements as boundary conditions can corrupt the material property estimates.[9–11] On the other hand, efficient direct methods to compute material properties from full vector displacement field measurements exist.[12–14] The present contribution uses an inverse problem approach that aims to improve lateral displacement estimates from ultrasound measurements.

Much research has been performed in order to recover better lateral displacements both in the UEI and the ultrasound blood flow imaging field. There are four main approaches that dominate the literature. In [15–25], the authors use various forms of interpolation within standard or modified block-matching algorithms to estimate subsample lateral displacements. These approaches try to make the best use of available data as is, but do not overcome the fundamental difficulty posed by the anisotropic PSF. An alternative approach is to use beam steering to measure displacements at arbitrary angles, and then use the arbitrary angled displacement estimates to determine the axial and lateral displacements as described in [26–31]. This method uses additional image data to overcome the limitations of a single PSF. A third approach, pursued by [32–35], is to use novel beam forming techniques to modulate the RF signal in the lateral direction so that it contains more phase information. This approach tries to reduce the anisotropy of the PSF. A fourth approach, pursued by [11,36–38], is to constrain the displacements to satisfy the incompressibility condition. Many different types of soft tissues can be assumed to be incompressible because they are mostly composed of water. A fifth approach, pursued in [39], is to constrain the displacements to satisfy the momentum equations.

Among the existing approaches, the incompressibility processing method introduced in [36], and the momentum equation based processing method [39] are arguably the most similar in spirit to the method presented here to reconstruct better lateral displacements. In [36], the authors assumed that the deformation is linear, incompressible and plane strain in character. With these assumptions, they formulated a minimization problem where the objective was to find lateral displacements that minimized the error between the measurements and predicted lateral displacements that satisfied the modelling assumptions. The incompressibility processing described in [38] is basically the same as the method developed in [36], with the addition of weighting functions to the lateral displacement reconstruction equations. These weighting functions are based on the absolute value and the gradient of the cross-correlation function. In [37], the method developed in [36] is extended to accommodate large plane strain deformations. Finally, in [11], the incompressibility constraint is used to improve the lateral and elevational displacements for a 3D deformation.

In [39], the authors constrain the measured displacements to satisfy the plane strain elasticity equations. Using only the measured axial displacements, and an $H^{\infty }$ filtering strategy, they are able to reconstruct the full displacement field. The drawback of this approach is that it requires a priori knowledge of the material parameters’ distribution.

In this paper, we also use the momentum equation to reconstruct the lateral displacements. The approach assumes only that the material property distribution is well-represented as piecewise constant. We will be using a spatial regularization term within an inverse problem formulation to adaptively enforce, and locally relax, a special form of the momentum conservation equation on the measured displacement field. Special weighting functions will be used to place more emphasis on the axial component of the displacement field.

The rest of the paper is arranged as follows: Section 2 describes the new formulation to estimate better lateral displacements; we refer to this method as ‘SPREME’, which stands for SParse RElaxation of the Momentum Equation. Section 3 describes the simulation experiments used to test and verify the displacement filtering algorithm. Section 4 describes the phantom experiments used to validate the performance of the displacement filtering algorithm on ultrasound measured data. Section 5 describes the in-vivo experiments used to further test the performance of the displacement filtering algorithm on patient data. Finally, a brief summary and conclusion is given in Section 6 along with possible future directions of the work. The appendix describes a simpler formulation that produces apparently improved displacement and strain fields. Despite appearances, however, the resulting strain and displacement fields are non-physical.

2. Formulation

We suppose we are given a measured displacements in a 2D region $\mathrm{\Omega }$. We further suppose the measurements of the axial displacement, $u_y$, are significantly more precise and accurate than the measurements of the lateral displacement. We make use of the assumption that within the observed plane, the body’s deformation may be well approximated by the plane stress approximation. These assumptions lead to the following problem statement:

Problem 2.1:

Given 2D measured displacement ${\mathit{u}}_m(x,y)$, for all $\mathit{x}\in \mathrm{\Omega }$, find the 2D displacement vector $\mathit{u}(x,y)$, and the 2D linearized strain tensor $\mathit{ϵ}(x,y)$ such that:(1) $$\begin{array}{c}\hfill {\mathit{\pi }}_o[\mathit{u}]=\frac{1}{2}||({\mathit{u}}_m-\mathit{u}){||}^2=\frac{1}{2}{\int }_{\mathrm{\Omega }}{({\mathit{u}}_m-\mathit{u})}^{⊺}·\mathit{T}({\mathit{u}}_m-\mathit{u})\mathrm{d}\mathrm{\Omega }\end{array}$$(1)

is minimized, and:(2) $$\begin{array}{c}\hfill \mathit{ϵ}={\mathrm{\nabla }}^s\mathit{u}\end{array}$$(2)

and:(3) $$\begin{array}{c}\hfill \mathrm{\nabla }·\mathbf{A}=0\text{a.e.}\text{in}\mathrm{\Omega }\end{array}$$(3)

where:(4) $$\begin{array}{c}\hfill \mathbf{A}(\mathit{ϵ})=2\text{tr}(\mathit{ϵ})\mathit{I}+2\mathit{ϵ}.\end{array}$$(4)

and:(5) $$\begin{array}{c}\hfill {\mathrm{\nabla }}^s\mathit{u}=\frac{1}{2}\left(\mathrm{\nabla }\mathit{u}+{(\mathrm{\nabla }\mathit{u})}^{badhbox}\right).\end{array}$$(5)

With diagonal $\mathit{T}$, Equation (1(1) $$\begin{array}{c}\hfill {\mathit{\pi }}_o[\mathit{u}]=\frac{1}{2}||({\mathit{u}}_m-\mathit{u}){||}^2=\frac{1}{2}{\int }_{\mathrm{\Omega }}{({\mathit{u}}_m-\mathit{u})}^{⊺}·\mathit{T}({\mathit{u}}_m-\mathit{u})\mathrm{d}\mathrm{\Omega }\end{array}$$(1) ) may be written explicitly as:(6) $$\begin{array}{c}\hfill {\mathit{\pi }}_o[\mathit{u}]=\frac{1}{2}{\int }_{\mathrm{\Omega }}[{\text{T}}_{xx}{(u_{(m)x}-u_x)}^2+{\text{T}}_{yy}{(u_{(m)y}-u_y)}^2]\mathrm{d}\mathrm{\Omega }.\end{array}$$(6)

In Equation (6(6) $$\begin{array}{c}\hfill {\mathit{\pi }}_o[\mathit{u}]=\frac{1}{2}{\int }_{\mathrm{\Omega }}[{\text{T}}_{xx}{(u_{(m)x}-u_x)}^2+{\text{T}}_{yy}{(u_{(m)y}-u_y)}^2]\mathrm{d}\mathrm{\Omega }.\end{array}$$(6) ), the weights (T${}_{xx}$, T${}_{yy}$) allow more importance to be placed on the accurate component of the displacement field while calculating the predicted displacement field. In Equation (3(3) $$\begin{array}{c}\hfill \mathrm{\nabla }·\mathbf{A}=0\text{a.e.}\text{in}\mathrm{\Omega }\end{array}$$(3) ), a.e. means almost everywhere, and implies that the constraint is enforced everywhere in the domain except on a set of measure zero (i.e. along curves and at isolated points).

To motivate constraint (3(3) $$\begin{array}{c}\hfill \mathrm{\nabla }·\mathbf{A}=0\text{a.e.}\text{in}\mathrm{\Omega }\end{array}$$(3) ), we consider a thin sheet in the (x, y) plane made of solid material that is linearly elastic, isotropic, incompressible and simply connected. We further assume that the deformation of the material is small, quasistatic and approximately plane stress in character. These modelling assumptions are justified because our application of interest is in the ultrasound elastography field, where the goal is to determine the linear elastic properties of soft tissue, typically modelled as incompressible because of their usually large water content, from small displacement measurements. The assumptions lead to the following definitions for the linear strain $\mathit{ϵ}$, and the Cauchy stress $\mathit{\sigma }$ [40,41]:(7) $$\begin{array}{cc}\hfill {\mathit{ϵ}}_{ij}& =\frac{1}{2}({\mathit{\partial }}_i{u}_j+{\mathit{\partial }}_j{u}_i),i,j=1,2\hfill \end{array}$$(7) (8) $$\begin{array}{cc}\hfill \mathit{\sigma }& =-p\mathit{I}+2\mathit{\mu }\mathit{ϵ},\hfill \end{array}$$(8)

where ${\mathit{\partial }}_i=\frac{\mathit{\partial }}{\mathit{\partial }{x}_i}$, $\mathit{\mu }(x,y)$ is the shear modulus distribution, and p(x, y) is the hydrostatic pressure field. For a plane stress deformation:(9) $$\begin{array}{c}\hfill {\mathit{\sigma }}_{zz}=-p+2\mathit{\mu }{\mathit{ϵ}}_{zz}=0.\end{array}$$(9)

For an incompressible material:(10) $$\begin{array}{c}\hfill {\mathit{ϵ}}_{xx}+{\mathit{ϵ}}_{yy}+{\mathit{ϵ}}_{zz}=0.\end{array}$$(10)

Substituting (10(10) $$\begin{array}{c}\hfill {\mathit{ϵ}}_{xx}+{\mathit{ϵ}}_{yy}+{\mathit{ϵ}}_{zz}=0.\end{array}$$(10) ) into (9(9) $$\begin{array}{c}\hfill {\mathit{\sigma }}_{zz}=-p+2\mathit{\mu }{\mathit{ϵ}}_{zz}=0.\end{array}$$(9) ) gives:(11) $$\begin{array}{c}\hfill p=-2\mathit{\mu }({\mathit{ϵ}}_{xx}+{\mathit{ϵ}}_{yy}).\end{array}$$(11)

Using (11(11) $$\begin{array}{c}\hfill p=-2\mathit{\mu }({\mathit{ϵ}}_{xx}+{\mathit{ϵ}}_{yy}).\end{array}$$(11) ) in (8(8) $$\begin{array}{cc}\hfill \mathit{\sigma }& =-p\mathit{I}+2\mathit{\mu }\mathit{ϵ},\hfill \end{array}$$(8) ) yields:(12) $$\begin{array}{c}\hfill \mathit{\sigma }=2\mathit{\mu }({\mathit{ϵ}}_{xx}+{\mathit{ϵ}}_{yy})\mathit{I}+2\mathit{\mu }\mathit{ϵ}=\mathit{\mu }\mathbf{A}.\end{array}$$(12)

The equilibrium equation, with no body force, is $\mathrm{\nabla }·\mathit{\sigma }=0$, or using (12(12) $$\begin{array}{c}\hfill \mathit{\sigma }=2\mathit{\mu }({\mathit{ϵ}}_{xx}+{\mathit{ϵ}}_{yy})\mathit{I}+2\mathit{\mu }\mathit{ϵ}=\mathit{\mu }\mathbf{A}.\end{array}$$(12) ):(13) $$\begin{array}{c}\hfill \mathrm{\nabla }·(\mathit{\mu }\mathbf{A})=0\text{in}\mathrm{\Omega }.\end{array}$$(13)

If $\mathit{\mu }(x,y)$ is piecewise constant, then $\mathrm{\nabla }·\mathbf{A}=0$ a.e. in $\mathrm{\Omega }$.

2.1. Uniqueness of solution

Given an axial displacement measurement $u_y$, the conditions necessary to guarantee the uniqueness of the lateral displacements computed using the equilibrium equation constraints will be derived in this section. This derivation begins with the assumption that the material has a piecewise homogeneous shear modulus distribution. This means that the material’s domain, $\mathrm{\Omega }$, can be divided into n different subdomains, ${\mathrm{\Omega }}^{(n)}$, and the shear modulus in each subdomain is constant. The shear modulus is not constant on the interface between the subdomains, however. An illustration of this type of domain, for a special case when $n=2$, is shown in Figure 1. The conditions necessary for uniqueness of the lateral displacements in each subdomain will be derived next.

2.1.1. Solution in a single subdomain

Lemma 2.1:

Given $u_y(x,y)$ in subdomain ${\mathrm{\Omega }}^{(n)}$, in which $\mathit{\mu }\equiv$ constant, then the general solution of Equations (2(2) $$\begin{array}{c}\hfill \mathit{ϵ}={\mathrm{\nabla }}^s\mathit{u}\end{array}$$(2) )–(5(5) $$\begin{array}{c}\hfill {\mathrm{\nabla }}^s\mathit{u}=\frac{1}{2}\left(\mathrm{\nabla }\mathit{u}+{(\mathrm{\nabla }\mathit{u})}^{badhbox}\right).\end{array}$$(5) ) for $u_x(x,y)$, with $(x,y)\in {\mathrm{\Omega }}^{(n)}$ is:(14) $$\begin{array}{c}\hfill u_x(x,y)=u_x^p(x,y)+c_{{}_1}(x^2-4y^2)+c_{{}_2}x+c_{{}_3}y+c_{{}_4}.\end{array}$$(14)

Proof:

Since $\mathit{\mu }=constant$ in ${\mathrm{\Omega }}^{(n)}$, we may expand Equation (3(3) $$\begin{array}{c}\hfill \mathrm{\nabla }·\mathbf{A}=0\text{a.e.}\text{in}\mathrm{\Omega }\end{array}$$(3) ) in terms of displacements via Equations (4(4) $$\begin{array}{c}\hfill \mathbf{A}(\mathit{ϵ})=2\text{tr}(\mathit{ϵ})\mathit{I}+2\mathit{ϵ}.\end{array}$$(4) ) and (2(2) $$\begin{array}{c}\hfill \mathit{ϵ}={\mathrm{\nabla }}^s\mathit{u}\end{array}$$(2) ), resulting in the following equations:(15) $$\begin{array}{cc}\hfill -4u_{x,xx}-u_{x,yy}& =3u_{y,yx},\hfill \end{array}$$(15) (16) $$\begin{array}{cc}\hfill -3u_{x,xy}& =u_{y,xx}+4u_{y,yy}.\hfill \end{array}$$(16)

All the terms on the right-hand side of both these equations are supposed to be known since by assumption here, $u_y$ is known precisely. The terms on the left-hand side involving the lateral displacements, $u_x$, are unknowns. We now show that these equations lead to a solution for $u_x$ that has at most four unknown constants. We begin by splitting the total solution into homogeneous and particular parts. To that end, we let:(17) $$\begin{array}{c}\hfill u_x=u_x^p+u_x^h.\end{array}$$(17)

Here, $u_x^p$ is any particular solution of (15(15) $$\begin{array}{cc}\hfill -4u_{x,xx}-u_{x,yy}& =3u_{y,yx},\hfill \end{array}$$(15) ) and (16(16) $$\begin{array}{cc}\hfill -3u_{x,xy}& =u_{y,xx}+4u_{y,yy}.\hfill \end{array}$$(16) ), which, by assumption, has no free constants or undetermined functions. If they appear in the process of solving (15(15) $$\begin{array}{cc}\hfill -4u_{x,xx}-u_{x,yy}& =3u_{y,yx},\hfill \end{array}$$(15) ) and (16(16) $$\begin{array}{cc}\hfill -3u_{x,xy}& =u_{y,xx}+4u_{y,yy}.\hfill \end{array}$$(16) ), then they shall be set to arbitrary values (e.g. zero) without loss of generality. Any part of the general solution (17(17) $$\begin{array}{c}\hfill u_x=u_x^p+u_x^h.\end{array}$$(17) ) that is undetermined by $u_y$ is thus contained in $u_x^h$, the homogeneous solution of (15(15) $$\begin{array}{cc}\hfill -4u_{x,xx}-u_{x,yy}& =3u_{y,yx},\hfill \end{array}$$(15) ) and (16(16) $$\begin{array}{cc}\hfill -3u_{x,xy}& =u_{y,xx}+4u_{y,yy}.\hfill \end{array}$$(16) ).

The homogeneous solution satisfies:(18) $$\begin{array}{cc}\hfill 4u_{x,xx}^h+u_{x,yy}^h& =0,\hfill \end{array}$$(18) (19) $$\begin{array}{cc}\hfill 3u_{x,xy}^h& =0.\hfill \end{array}$$(19)

Integrating (19(19) $$\begin{array}{cc}\hfill 3u_{x,xy}^h& =0.\hfill \end{array}$$(19) ) twice yields:(20) $$\begin{array}{c}\hfill u_x^h=f(x)+g(y).\end{array}$$(20)

Substituting (20(20) $$\begin{array}{c}\hfill u_x^h=f(x)+g(y).\end{array}$$(20) ) into (18(18) $$\begin{array}{cc}\hfill 4u_{x,xx}^h+u_{x,yy}^h& =0,\hfill \end{array}$$(18) ) gives:(21) $$\begin{array}{c}\hfill 4f^{″}(x)+g^{″}(y)=0.\end{array}$$(21)

Equation (21(21) $$\begin{array}{c}\hfill 4f^{″}(x)+g^{″}(y)=0.\end{array}$$(21) ) implies:(22) $$\begin{array}{c}\hfill 4f^{″}(x)=-g^{″}(y)=8c_{{}_1}.\end{array}$$(22)

Integrating (22(22) $$\begin{array}{c}\hfill 4f^{″}(x)=-g^{″}(y)=8c_{{}_1}.\end{array}$$(22) ) for f(x) and g(y) yields:(23) $$\begin{array}{cc}\hfill f(x)& =c_{{}_1}{x}^2+c_{{}_2}x+c_{{}_4}^{\prime },\hfill \end{array}$$(23) (24) $$\begin{array}{cc}\hfill g(y)& =-4c_{{}_1}{y}^2+c_{{}_3}y+c_{{}_5}^{\prime }.\hfill \end{array}$$(24)

Combining (23(23) $$\begin{array}{cc}\hfill f(x)& =c_{{}_1}{x}^2+c_{{}_2}x+c_{{}_4}^{\prime },\hfill \end{array}$$(23) ), (24(24) $$\begin{array}{cc}\hfill g(y)& =-4c_{{}_1}{y}^2+c_{{}_3}y+c_{{}_5}^{\prime }.\hfill \end{array}$$(24) ), in (20(20) $$\begin{array}{c}\hfill u_x^h=f(x)+g(y).\end{array}$$(20) ) gives:(25) $$\begin{array}{c}\hfill u_x^h(x,y)=c_{{}_1}(x^2-4y^2)+c_{{}_2}x+c_{{}_3}y+c_{{}_4},\end{array}$$(25)

where $c_{{}_4}=c_{{}_4}^{\prime }+c_{{}_5}^{\prime }$. Therefore, the total solution is of the form (14(14) $$\begin{array}{c}\hfill u_x(x,y)=u_x^p(x,y)+c_{{}_1}(x^2-4y^2)+c_{{}_2}x+c_{{}_3}y+c_{{}_4}.\end{array}$$(14) ). $\square$

Equation (14(14) $$\begin{array}{c}\hfill u_x(x,y)=u_x^p(x,y)+c_{{}_1}(x^2-4y^2)+c_{{}_2}x+c_{{}_3}y+c_{{}_4}.\end{array}$$(14) ) shows that within any homogeneous subregion of the body, given the axial displacement $u_y(x,y)$, the lateral displacement $u_x(x,y)$ has at most four undetermined constants.

2.1.2. Solution in connected subdomains covering entire domain of interest

Interestingly, the form of the homogeneous solution is completely independent of the given $u_y(x,y)$. Thus, we may conclude almost immediately:

Lemma 2.2:

Given $u_y(x,y)$ in $\mathrm{\Omega }=\cup {\mathrm{\Omega }}^{(n)}$, in which $\mathit{\mu }(x,y)$ is piecewise constant, such that(26) $$\begin{array}{c}\hfill \mathit{\mu }(x,y)={\mathit{\mu }}^{(n)}\text{if}(x,y)\in {\mathrm{\Omega }}^{(n)},\end{array}$$(26)

then the general solution of Equations (2(2) $$\begin{array}{c}\hfill \mathit{ϵ}={\mathrm{\nabla }}^s\mathit{u}\end{array}$$(2) )–(5(5) $$\begin{array}{c}\hfill {\mathrm{\nabla }}^s\mathit{u}=\frac{1}{2}\left(\mathrm{\nabla }\mathit{u}+{(\mathrm{\nabla }\mathit{u})}^{badhbox}\right).\end{array}$$(5) ) for $u_x(x,y)$ with $(x,y)\in \mathrm{\Omega }=\cup {\mathrm{\Omega }}^{(n)}$ is:(27) $$\begin{array}{c}\hfill u_x(x,y)=u_x^p(x,y)+c_{{}_1}(x^2-4y^2)+c_{{}_2}x+c_{{}_3}y+c_{{}_4}.\end{array}$$(27)

Here, $u_x^p(x,y)$ is determined entirely by $u_y(x,y)$, while $c_{{}_1},\dots ,c_{{}_4}$ are undetermined constants.

Proof:

We consider several contiguous subregions, each of which is homogeneous. From Lemma 2.1, the solution within subregion n is given by:(28) $$\begin{array}{c}\hfill u_x^{(n)}(x,y)=u_x^{(n)p}(x,y)+c_{{}_1}^{(n)}(x^2-4y^2)+c_{{}_2}^{(n)}x+c_{{}_3}^{(n)}y+c_{{}_4}^{(n)}.\end{array}$$(28)

If we consider the displacement $u_x(x,y)$ to be continuous at a non-special interface¹ between subregion n and adjacent subregion m, and we assume that $u_x^{(n)p}=u_x^{(m)p}$ is continuous on the interface, then we conclude:(29) $$\begin{array}{c}\hfill c_{{}_j}^{(n)}=c_{{}_j}^{(m)}.\end{array}$$(29)

To see this, note that $u_x^{(n)}=u_x^{(m)}$, and $u_x^{(n)p}=u_x^{(m)p}$ leads to the equation:(30) $$\begin{array}{c}\hfill \left[\begin{array}{cccc}{x}^2-4y^2& x& y& 1\end{array}\right]\left\{\begin{array}{c}{c}_{{}_1}^{\ast }\\ c_{{}_2}^{\ast }\\ c_{{}_3}^{\ast }\\ c_{{}_4}^{\ast }\end{array}\right\}=\left\{\begin{array}{c}0\end{array}\right\}\end{array}$$(30)

with $c_{{}_j}^{\ast }=c_{{}_j}^{(n)}-c_{{}_j}^{(m)}$. We note that (30(30) $$\begin{array}{c}\hfill \left[\begin{array}{cccc}{x}^2-4y^2& x& y& 1\end{array}\right]\left\{\begin{array}{c}{c}_{{}_1}^{\ast }\\ c_{{}_2}^{\ast }\\ c_{{}_3}^{\ast }\\ c_{{}_4}^{\ast }\end{array}\right\}=\left\{\begin{array}{c}0\end{array}\right\}\end{array}$$(30) ) holds at every point on the interface. Provided that the functions $x^2-4y^2$, x, y, 1 are linearly independent on the interfaces, then the only solution of (30(30) $$\begin{array}{c}\hfill \left[\begin{array}{cccc}{x}^2-4y^2& x& y& 1\end{array}\right]\left\{\begin{array}{c}{c}_{{}_1}^{\ast }\\ c_{{}_2}^{\ast }\\ c_{{}_3}^{\ast }\\ c_{{}_4}^{\ast }\end{array}\right\}=\left\{\begin{array}{c}0\end{array}\right\}\end{array}$$(30) ) is $c_{{}_j}^{\ast }=0$. Interfaces on which these functions are linearly dependent are special cases, and we call them ‘special interfaces’.

Since the $c_j$’s take the same value in every subregion, we may consider the homogeneous functions to be defined globally. Thus, the global solution is:(31) $$\begin{array}{c}\hfill u_x(x,y)=u_x^p(x,y)+c_{{}_1}(x^2-4y^2)+c_{{}_2}x+c_{{}_3}y+c_{{}_4}.\end{array}$$(31) $\square$

Lemma 2.2 and Equation (31(31) $$\begin{array}{c}\hfill u_x(x,y)=u_x^p(x,y)+c_{{}_1}(x^2-4y^2)+c_{{}_2}x+c_{{}_3}y+c_{{}_4}.\end{array}$$(31) ) imply that by knowing the axial displacement throughout the domain in a piecewise homogeneous material, the lateral displacement may be determined up to four coefficients. Hence, if sufficient information exists in the measurements to determine these four coefficients accurately, then the lateral displacement field may be determined accurately throughout the domain.

In the next section, we derive a variational formulation that permits this estimation in a natural way.

2.2. SPREME variational formulation

In practice, we shall enforce constraint (3(3) $$\begin{array}{c}\hfill \mathrm{\nabla }·\mathbf{A}=0\text{a.e.}\text{in}\mathrm{\Omega }\end{array}$$(3) ) approximately as a penalty, in the form of a regularization term. To determine the appropriate form of that term, we refer to Figure 1 and consider a continuous $\mathit{\mu }$ of the form:(32) $$\begin{array}{c}\hfill \mathit{\mu }=\left\{\begin{array}{cc}{\mathit{\mu }}_1\hfill & x\in {\mathrm{\Omega }}_1\hfill \\ {\mathit{\mu }}_2\hfill & x\in {\mathrm{\Omega }}_2\hfill \\ \text{continuous monotone}\hfill & x\in {\mathrm{\Omega }}_c\hfill \end{array}\right.\end{array}$$(32)

Figure 1. A piecewise continuous distribution of elastic modulus represented as the limit of a continuous distribution. The function $\mathit{\mu }(x,y)$ has the constant value ${\mathit{\mu }}_1$ for $(x,y)\in {\mathrm{\Omega }}_1$, has the constant value ${\mathit{\mu }}_2$ for $(x,y)\in {\mathrm{\Omega }}_2$ and is continuous and monotone in a narrow intermediate region, ${\mathrm{\Omega }}_c$.

where ${\mathrm{\Omega }}_c$ is the set of all points within $\frac{h}{2}$ from curve C. We assume there is a constant $C_1$ so that $|\mathrm{\nabla }\mathit{\mu }|\le \frac{C_1}{h}$ for all $x\in {\mathrm{\Omega }}_c.$

Suppose $\mathrm{\nabla }·(\mathit{\mu }\mathbf{A})=0$ everywhere in $\mathrm{\Omega }$, then:(33) $$\begin{array}{c}\hfill \mathrm{\nabla }·(\mathit{\mu }\mathbf{A})=\left\{\begin{array}{cc}{\mathit{\mu }}_1\mathrm{\nabla }·\mathbf{A}=0\hfill & x\in {\mathrm{\Omega }}_1\hfill \\ {\mathit{\mu }}_2\mathrm{\nabla }·\mathbf{A}=0\hfill & x\in {\mathrm{\Omega }}_2\hfill \\ \mathit{\mu }\mathrm{\nabla }·\mathbf{A}+\mathbf{A}\mathrm{\nabla }\mathit{\mu }=0\hfill & x\in {\mathrm{\Omega }}_c\hfill \end{array}\right..\end{array}$$(33)

In order to form the penalty term, we assume (33(33) $$\begin{array}{c}\hfill \mathrm{\nabla }·(\mathit{\mu }\mathbf{A})=\left\{\begin{array}{cc}{\mathit{\mu }}_1\mathrm{\nabla }·\mathbf{A}=0\hfill & x\in {\mathrm{\Omega }}_1\hfill \\ {\mathit{\mu }}_2\mathrm{\nabla }·\mathbf{A}=0\hfill & x\in {\mathrm{\Omega }}_2\hfill \\ \mathit{\mu }\mathrm{\nabla }·\mathbf{A}+\mathbf{A}\mathrm{\nabla }\mathit{\mu }=0\hfill & x\in {\mathrm{\Omega }}_c\hfill \end{array}\right..\end{array}$$(33) ) is satisfied and consider:(34) $$\begin{array}{cc}\hfill {\int }_{\mathrm{\Omega }}{|\mathrm{\nabla }·\mathbf{A}|}^{\mathit{\alpha }}\mathrm{d}\mathrm{\Omega }& ={\int }_{{\mathrm{\Omega }}_c}{\left|\mathbf{A}\frac{\mathrm{\nabla }\mathit{\mu }}{\mathit{\mu }}\right|}^{\mathit{\alpha }}\mathrm{d}{\mathrm{\Omega }}_c\le {\int }_{{\mathrm{\Omega }}_c}\frac{{|\mathbf{A}|}^{\mathit{\alpha }}{|\mathrm{\nabla }\mathit{\mu }|}^{\mathit{\alpha }}}{{({\mathit{\mu }}_{\min })}^{\mathit{\alpha }}}\mathrm{d}{\mathrm{\Omega }}_c\le {\int }_{{\mathrm{\Omega }}_c}{\left(\frac{C_2}{h}\right)}^{\mathit{\alpha }}\mathrm{d}{\mathrm{\Omega }}_c\hfill \\ \hfill & \approx \hfill & \hfill h{L}_{(c)}{\left(\frac{C_2}{h}\right)}^{\mathit{\alpha }}=L_{(c)}{C}_2^{\mathit{\alpha }}{h}^{1-\mathit{\alpha }}\end{array}$$(34)

where ${\mathit{\mu }}_{min}=$ min$\{{\mathit{\mu }}_1,{\mathit{\mu }}_2\}$ is the minimum value of $\mathit{\mu }$ in ${\mathrm{\Omega }}_c$, $C_2=(C_1|\mathbf{A}|)/{\mathit{\mu }}_{min}$, $|\mathbf{A}|$ is the largest eigenvalue of $\mathbf{A}$ and $L_{(c)}$ is the perimeter of curve C; see Figure 1.

We see in (34(34) $$\begin{array}{cc}\hfill {\int }_{\mathrm{\Omega }}{|\mathrm{\nabla }·\mathbf{A}|}^{\mathit{\alpha }}\mathrm{d}\mathrm{\Omega }& ={\int }_{{\mathrm{\Omega }}_c}{\left|\mathbf{A}\frac{\mathrm{\nabla }\mathit{\mu }}{\mathit{\mu }}\right|}^{\mathit{\alpha }}\mathrm{d}{\mathrm{\Omega }}_c\le {\int }_{{\mathrm{\Omega }}_c}\frac{{|\mathbf{A}|}^{\mathit{\alpha }}{|\mathrm{\nabla }\mathit{\mu }|}^{\mathit{\alpha }}}{{({\mathit{\mu }}_{\min })}^{\mathit{\alpha }}}\mathrm{d}{\mathrm{\Omega }}_c\le {\int }_{{\mathrm{\Omega }}_c}{\left(\frac{C_2}{h}\right)}^{\mathit{\alpha }}\mathrm{d}{\mathrm{\Omega }}_c\hfill \\ \hfill & \approx \hfill & \hfill h{L}_{(c)}{\left(\frac{C_2}{h}\right)}^{\mathit{\alpha }}=L_{(c)}{C}_2^{\mathit{\alpha }}{h}^{1-\mathit{\alpha }}\end{array}$$(34) ) that for piecewise constant $\mathit{\mu }$ (i.e. when $h\to 0$),(35) $$\begin{array}{c}\hfill {\int }_{\mathrm{\Omega }}{|\mathrm{\nabla }·\mathbf{A}|}^{\mathit{\alpha }}\mathrm{d}\mathrm{\Omega }=0\text{for}0<\mathit{\alpha }<1.\end{array}$$(35)

Equation (35(35) $$\begin{array}{c}\hfill {\int }_{\mathrm{\Omega }}{|\mathrm{\nabla }·\mathbf{A}|}^{\mathit{\alpha }}\mathrm{d}\mathrm{\Omega }=0\text{for}0<\mathit{\alpha }<1.\end{array}$$(35) ) shows that we may enforce the constraint $\mathrm{\nabla }·\mathbf{A}=0a.e.$ by choosing to enforce the integral condition (35(35) $$\begin{array}{c}\hfill {\int }_{\mathrm{\Omega }}{|\mathrm{\nabla }·\mathbf{A}|}^{\mathit{\alpha }}\mathrm{d}\mathrm{\Omega }=0\text{for}0<\mathit{\alpha }<1.\end{array}$$(35) ) with $0<\mathit{\alpha }<1$.

The derivative of the absolute value function is undefined at the origin, so in practice, we introduce a small positive constant, $\mathit{\delta }$, and approximate (35(35) $$\begin{array}{c}\hfill {\int }_{\mathrm{\Omega }}{|\mathrm{\nabla }·\mathbf{A}|}^{\mathit{\alpha }}\mathrm{d}\mathrm{\Omega }=0\text{for}0<\mathit{\alpha }<1.\end{array}$$(35) ) as:(36) $$\begin{array}{c}\hfill {\int }_{\mathrm{\Omega }}{|\mathrm{\nabla }·\mathbf{A}(\mathit{ϵ})|}^{\mathit{\alpha }}\mathrm{d}\mathrm{\Omega }\approx {\int }_{\mathrm{\Omega }}\frac{{\left(\mathrm{\nabla }·\mathbf{A}({\mathit{ϵ}}^k)\right)}^2}{[{\left(\mathrm{\nabla }·\mathbf{A}({\mathit{ϵ}}^{{}^{(k-1)}})\right)}^2+\mathit{\delta }{]}^n}\mathrm{d}\mathrm{\Omega }\end{array}$$(36)

where:(37) $$\begin{array}{c}\hfill n=1-\frac{\mathit{\alpha }}{2}.\end{array}$$(37)

Note that $0<\mathit{\alpha }<1⟺1>n>\frac{1}{2}$.

In Equation (36(36) $$\begin{array}{c}\hfill {\int }_{\mathrm{\Omega }}{|\mathrm{\nabla }·\mathbf{A}(\mathit{ϵ})|}^{\mathit{\alpha }}\mathrm{d}\mathrm{\Omega }\approx {\int }_{\mathrm{\Omega }}\frac{{\left(\mathrm{\nabla }·\mathbf{A}({\mathit{ϵ}}^k)\right)}^2}{[{\left(\mathrm{\nabla }·\mathbf{A}({\mathit{ϵ}}^{{}^{(k-1)}})\right)}^2+\mathit{\delta }{]}^n}\mathrm{d}\mathrm{\Omega }\end{array}$$(36) ), k denotes an iteration counter. As k approaches infinity, we expect ${\mathit{ϵ}}^{{}^{(k-1)}}$ to approach ${\mathit{ϵ}}^k$. In this limit, for $\mathit{\delta }=0$, the approximation in (36(36) $$\begin{array}{c}\hfill {\int }_{\mathrm{\Omega }}{|\mathrm{\nabla }·\mathbf{A}(\mathit{ϵ})|}^{\mathit{\alpha }}\mathrm{d}\mathrm{\Omega }\approx {\int }_{\mathrm{\Omega }}\frac{{\left(\mathrm{\nabla }·\mathbf{A}({\mathit{ϵ}}^k)\right)}^2}{[{\left(\mathrm{\nabla }·\mathbf{A}({\mathit{ϵ}}^{{}^{(k-1)}})\right)}^2+\mathit{\delta }{]}^n}\mathrm{d}\mathrm{\Omega }\end{array}$$(36) ) is then exact. The $\mathit{\delta }$ in the denominator is a small numerical constant that prevents the denominator from being zero when $\mathrm{\nabla }·\mathbf{A}=0.$

Using (36(36) $$\begin{array}{c}\hfill {\int }_{\mathrm{\Omega }}{|\mathrm{\nabla }·\mathbf{A}(\mathit{ϵ})|}^{\mathit{\alpha }}\mathrm{d}\mathrm{\Omega }\approx {\int }_{\mathrm{\Omega }}\frac{{\left(\mathrm{\nabla }·\mathbf{A}({\mathit{ϵ}}^k)\right)}^2}{[{\left(\mathrm{\nabla }·\mathbf{A}({\mathit{ϵ}}^{{}^{(k-1)}})\right)}^2+\mathit{\delta }{]}^n}\mathrm{d}\mathrm{\Omega }\end{array}$$(36) ) instead of (3(3) $$\begin{array}{c}\hfill \mathrm{\nabla }·\mathbf{A}=0\text{a.e.}\text{in}\mathrm{\Omega }\end{array}$$(3) ), we reformulate the problem as:

Problem 2.2:

Given ${\mathit{u}}_m(x,y)$, $\mathit{x}\in \mathrm{\Omega }$ and ${\mathit{ϵ}}^{(k-1)}(x,y)$, find ${\mathit{u}}^k$ and ${\mathit{ϵ}}^k$ that minimize:(38) $$\begin{array}{c}\hfill \mathit{\pi }[{\mathit{u}}^k,{\mathit{ϵ}}^k]={\mathit{\pi }}_{{}_O}+{\mathit{\pi }}_{{}_C}+{\mathit{\pi }}_{{}_R}\end{array}$$(38)

where:(39) $$\begin{array}{c}\hfill {\mathit{\pi }}_{{}_O}=\frac{1}{2}{\int }_{\mathrm{\Omega }}({\mathit{u}}_m-{\mathit{u}}^k)·\mathit{T}({\mathit{u}}_m-{\mathit{u}}^k)\mathrm{d}\mathrm{\Omega }\end{array}$$(39)

and:(40) $$\begin{array}{c}\hfill {\mathit{\pi }}_{{}_C}=\frac{\mathit{\beta }}{2}{\int }_{\mathrm{\Omega }}{({\mathit{ϵ}}^k-{\mathrm{\nabla }}^s{\mathit{u}}^k)}^2\mathrm{d}\mathrm{\Omega }\end{array}$$(40)

and:(41) $$\begin{array}{c}\hfill {\mathit{\pi }}_{{}_R}=\frac{1}{2}{\int }_{\mathrm{\Omega }}\mathit{\alpha }(\mathit{x}){\left(\mathrm{\nabla }·\mathbf{A}({\mathit{ϵ}}^k)\right)}^2\mathrm{d}\mathrm{\Omega }\end{array}$$(41)

and:(42) $$\begin{array}{c}\hfill \mathit{\alpha }(\mathit{x})=\frac{{\mathit{\alpha }}_o}{[{\left(\mathrm{\nabla }·\mathbf{A}({\mathit{ϵ}}^{{}^{(k-1)}})\right)}^2+\mathit{\delta }{]}^n}.\end{array}$$(42)

This yields a quadratic minimization problem at each iteration for $u^k$ and ${\mathit{ϵ}}^k$ (Note that the shear modulus $\mathit{\mu }$ is not required in the formulation). We initialize iterations using ${\mathit{ϵ}}^{(0)}=0$. Equation (38(38) $$\begin{array}{c}\hfill \mathit{\pi }[{\mathit{u}}^k,{\mathit{ϵ}}^k]={\mathit{\pi }}_{{}_O}+{\mathit{\pi }}_{{}_C}+{\mathit{\pi }}_{{}_R}\end{array}$$(38) ) contains a displacement data matching term (${\mathit{\pi }}_{{}_O}$), a compatibility term (${\mathit{\pi }}_{{}_C}$) and a regularization term (${\mathit{\pi }}_{{}_R}$). ${\mathit{\pi }}_{{}_O}$ minimizes the mismatch between the measured displacement field and the filtered displacement field and ${\mathit{\pi }}_{{}_C}$ constrains the calculated strains to satisfy the compatibility equation. $\mathit{\beta }$ is a constant that is used to regulate how strongly the compatibility condition is enforced. ${\mathit{\pi }}_{{}_R}$ is used to constrain the strains to satisfy the equilibrium equations almost everywhere in the domain. ${\mathit{\alpha }}_o$ is a constant used to regulate how strongly the equilibrium is enforced globally, and $\mathit{\alpha }(\mathit{x})$ may be interpreted as a spatially varying regularization constant that regulates how strongly the equilibrium equation is enforced locally.

The way the regularization, ${\mathit{\pi }}_{{}_R}$, works can be best explained with Figure 1. In this figure, where $\mathit{\mu }$ is constant (i.e. in the background (${\mathit{\mu }}_2$) and in the inclusion (${\mathit{\mu }}_1$)), $\mathrm{\nabla }·\mathbf{A}({\mathit{ϵ}}^{{}^{\left(k-1\right)}})$ will be small, meaning that $\mathit{\alpha }(\mathit{x})$ will be large, so the equilibrium equation will be strongly enforced in these regions. Therefore, $u_x$ can be confidently determined there up to the homogeneous solution. Where $\mathit{\mu }$ is not constant (i.e. in ${\mathrm{\Omega }}_c$, a thin region along the interface between the inclusion and the background), $\mathrm{\nabla }·\mathbf{A}({\mathit{ϵ}}^{{}^{\left(k-1\right)}})$ will be large, meaning that $\mathit{\alpha }(\mathit{x})$ will be small, so the equilibrium equation will not be strongly enforced. Thus, the strains are permitted to change rapidly there. Continuity of the lateral displacements is enforced there through the ${\mathit{\pi }}_{{}_C}$ term. This ‘spatially adaptive’ regularization thus allows us to enforce the equilibrium equation everywhere in the domain, except near material interfaces.

We note that if the goal is to evaluate the strain tensor only, then it is tempting to reformulate the problem without the displacement variable, $\mathit{u}$. Such a formulation is presented in Appendix 1. There, it is shown that the strains obtained from the formulation are not compatible, and hence not physical, for some simulation experiments.

Equation (38(38) $$\begin{array}{c}\hfill \mathit{\pi }[{\mathit{u}}^k,{\mathit{ϵ}}^k]={\mathit{\pi }}_{{}_O}+{\mathit{\pi }}_{{}_C}+{\mathit{\pi }}_{{}_R}\end{array}$$(38) ) is the SPREME functional. Our implementation minimizes $\mathit{\pi }[\mathit{u},\mathit{ϵ}]$ by the Galerkin finite element method (FEM) discretization. The weak form is derived first by computing the following directional derivatives:(43) $$\begin{array}{cc}\hfill D_{{}_{\mathit{u}}}\mathit{\pi }:\mathit{v}& =\frac{\mathrm{d}}{\mathrm{d}c}{|}_{c=0}\mathit{\pi }[\mathit{u}+c\mathit{v},\mathit{ϵ}]=0,\hfill \end{array}$$(43) (44) $$\begin{array}{cc}\hfill D_{{}_{\mathit{ϵ}}}\mathit{\pi }:\mathit{w}& =\frac{\mathrm{d}}{\mathrm{d}c}{|}_{c=0}\mathit{\pi }[\mathit{u},\mathit{ϵ}+c\mathit{w}]=0,\hfill \end{array}$$(44)

where $\mathit{v}$ and $\mathit{w}$ are admissible variations of $\mathit{u}$ and $\mathit{ϵ}$, respectively, and $(\mathit{u},\mathit{ϵ})=({\mathit{u}}^k,{\mathit{ϵ}}^k)$. We drop the k’s for the rest of the presentation to simplify notation. Taking the directional derivatives yields the following system of equations:(45) $$\begin{array}{cc}& {\int }_{\mathrm{\Omega }}\mathit{T}(\mathit{u}-{\mathit{u}}_m)\mathit{v}\mathrm{d}\mathrm{\Omega }-\mathit{\beta }{\int }_{\mathrm{\Omega }}\left[\mathit{ϵ}-{\mathrm{\nabla }}^s\mathit{u}\right](\mathrm{\nabla }\mathit{v})\mathrm{d}\mathrm{\Omega }=0\text{for all}\mathit{v},\hfill \end{array}$$(45) (46) $$\begin{array}{cc}& \mathit{\beta }{\int }_{\mathrm{\Omega }}\left[\mathit{ϵ}-{\mathrm{\nabla }}^s\mathit{u}\right]\mathit{w}\mathrm{d}\mathrm{\Omega }+{\int }_{\mathrm{\Omega }}\mathit{\alpha }(\mathit{x})\{\mathrm{\nabla }·\mathbf{A}(\mathit{w})\}·\{\mathrm{\nabla }·\mathbf{A}(\mathit{ϵ})\}\mathrm{d}\mathrm{\Omega }=0\text{for all}\mathit{w}.\hfill \end{array}$$(46)

Equations (45(45) $$\begin{array}{cc}& {\int }_{\mathrm{\Omega }}\mathit{T}(\mathit{u}-{\mathit{u}}_m)\mathit{v}\mathrm{d}\mathrm{\Omega }-\mathit{\beta }{\int }_{\mathrm{\Omega }}\left[\mathit{ϵ}-{\mathrm{\nabla }}^s\mathit{u}\right](\mathrm{\nabla }\mathit{v})\mathrm{d}\mathrm{\Omega }=0\text{for all}\mathit{v},\hfill \end{array}$$(45) ) and (46(46) $$\begin{array}{cc}& \mathit{\beta }{\int }_{\mathrm{\Omega }}\left[\mathit{ϵ}-{\mathrm{\nabla }}^s\mathit{u}\right]\mathit{w}\mathrm{d}\mathrm{\Omega }+{\int }_{\mathrm{\Omega }}\mathit{\alpha }(\mathit{x})\{\mathrm{\nabla }·\mathbf{A}(\mathit{w})\}·\{\mathrm{\nabla }·\mathbf{A}(\mathit{ϵ})\}\mathrm{d}\mathrm{\Omega }=0\text{for all}\mathit{w}.\hfill \end{array}$$(46) ) are discretized following standard FE methods, (e.g. [41]) with bilinear shape functions $N_a(\mathit{x})$ defined on quadrilateral elements. The shape functions have the special property that $N_a=1$ at node a and $N_a=0$ at all other nodes. This gives the following discrete approximations of the functions on each element:(47) $$\begin{array}{c}\hfill \mathit{v}\approx \sum _{a=1}^4{\mathit{v}}^a{N}_a(\mathit{x}),\mathit{u}\approx \sum _{b=1}^4{\mathit{u}}^b{N}_b(\mathit{x}),\mathit{w}\approx \sum _{a=1}^4{\mathit{w}}^a{N}_a(\mathit{x}),\mathit{ϵ}\approx \sum _{b=1}^4{\mathit{ϵ}}^b{N}_b(\mathit{x}).\end{array}$$(47)

Here, ${\mathit{u}}^b$, ${\mathit{ϵ}}^b$ represent the values of $\mathit{u}$ and $\mathit{ϵ}$ at node b.

This discretization of the variational equations results in an algebraic equation of the form $\mathit{K}\mathit{d}=\mathit{f}$, in which $\mathit{d}$ is a vector containing the unknown nodal values of $\mathit{u}$ and $\mathit{ϵ}$. A new finite element within an in-house FEM code was created to build and solve the system of linear equations. See [42] for more details on the finite element discretization.

3. Verification of computational implementation

The purpose of verification is to confirm that, given data that satisfy the modelling assumptions, a computational implementation produces the correct solution that is consistent with the assumed model. In this section, therefore, the results of two different simulation experiments conducted to test the implementation of the SPREME formulation will be described. For the two experiments performed, reference analytical displacement and strain fields are available, and they were used to evaluate the quality of reconstructed displacements and strains from the SPREME formulation.

3.1. Reference analytical solution

An exact analytical solution of the elasticity equations is used to generate reference displacement data. The solution corresponds to that of an infinite elastic sheet with a perfectly bonded circular inclusion. The sheet is subject to uniform stress at infinity. The problem was non-dimensionalized using the region of interest size as the reference length, and the background shear modulus as the reference stress. The shear modulus of the inclusion and background was ${\mathit{\mu }}_{{}_I}=3$, and ${\mathit{\mu }}_{{}_B}=1$, respectively. Poisson’s ratio was taken to be $\mathit{\nu }=0.5$, consistent with the incompressibility assumption. The radius of the inclusion was $a=0.25$. The solution to the elasticity equation was obtained using the methods described in [43].

3.2. Uniaxial tension test

Figure 2. Uniaxial and biaxial tension tests: The exact analytical solution for far-field tension applied to an infinite sheet with a bonded circular inclusion is used for verification of the computational implementation. In the figure, ROI indicates the region of interest over which displacement fields are assumed to have been ‘measured’.

Figure 3. Uniaxial tension test. Upper row: Exact noiseless displacement field for uniaxial tension experiment. These represent the target displacement fields. Lower row: Noise-corrupted displacement fields used for verification. The noise level in the lateral displacements is similar to that seen in physically measured displacement fields shown in Figure 10.

Figure 4. Uniaxial tension test results. Top row: Target noiseless strain fields from uniaxial tension test. Middle row: Strain fields calculated from noise-corrupted input displacement fields. Bottom row: Strain fields recovered from SPREME.

For the first experiment, the sheet was subject to uniaxial tension in the y direction as depicted in Figure 2. The stress at infinity was taken to be ${\mathit{\sigma }}_{yy}={\mathit{\sigma }}^{\infty }=6\times {10}^{-2}$. The displacement field corresponding to this deformation was evaluated in a subregion around the inclusion. This subregion was discretized into $50\times 50$ square elements in the x and y directions, and the axial and lateral displacements were evaluated at each of the 2601 nodes. This reference displacement field is shown in Figure 3. This target displacement field was differentiated by finite differences to calculate the target strains shown in Figure 4. Since the analytical expressions for the displacements are available, the target strains can be differentiated exactly rather than using finite differences. There is no analytical expression for the noise-corrupted displacement fields, however, and therefore we have chosen to treat both noisy and reference fields in the identical manner. These target fields will be used as a reference to benchmark the accuracy of the fields reconstructed with SPREME.

Our problem assumptions include that we are given accurate axial displacements and noisy lateral displacements. Therefore, Gaussian noise was added to the lateral displacement ($u_x$) in order to simulate ultrasound displacement measurements, and the axial displacement ($u_y$) was left unchanged. The corrupted displacement field is shown in Figure 3.

The corrupted displacement field was differentiated by finite differences to obtain the strains displayed in Figure 4. The ${\mathit{ϵ}}_{yy}$ strain component was unchanged. The features in the other strain components, however, were completely obscured by noise.

The corrupted displacement field was processed with SPREME using the following parameter values: ${\mathit{\alpha }}_o={10}^{-5}$, $\mathit{\beta }=1$, $\mathit{\delta }={10}^{-8}$, $n=0.5$, $T_{xx}={10}^{-9}$, $T_{yy}={10}^4$. These parameters were chosen by trial and error from extensive computational testing. This same trial and error method was used to choose the best parameters for all the other computational experiments performed. Six iterations were performed to allow the strains to converge. The strains from SPREME are displayed in the bottom row of Figure 4. The reconstructed strains are very smooth because of the regularization on the momentum equation. This regularization term is very similar to a total variation (TV) regularization which tends to penalize reconstructions with high oscillations, but allows reconstructions to have jumps that need to be present to fit the data.[44]

Accuracy of reconstructed strains and displacements may be seen qualitatively in line plots passing through the centre of the inclusion region, Figure 5. These plots show excellent agreement between the target strains and the reconstructed strains (referred to as SPREME Strain in the legend) in every component of the strain. We recall that the SPREME formulation has independent variables for both displacements and strains, and therefore we also use this figure to evaluate consistency of the SPREME displacements and strains. To see this, we use a finite difference approximation of the derivative of the reconstructed displacements to obtain another estimate of computed strains; these are labelled SPREME Disp in the legend, and are also in excellent agreement with the other quantities.

Figure 5. Uniaxial tension: Line plots through the axial, lateral and shear strain fields. The SPREME method produces independent estimates of displacement and strain. Thus, each graph shows three curves. One is the target strain. The second is the strain output from SPREME  labelled ‘SPREME Strain’. The third is the derivative of the output SPREME displacement, labelled ‘SPREME Disp’. Comparing the latter two provides a test for consistency between these fields.

Overall quantitative accuracy of the displacements (strains) computed with SPREME were evaluated for both individual components and for the total vector (tensor). The component error is defined as:(48) $$\begin{array}{c}\hfill \text{Strain error}=\frac{||{\mathit{ϵ}}_{ij}^t-{\mathit{ϵ}}_{ij}||}{||{\mathit{ϵ}}_{ij}^t||}\approx \frac{\sqrt{\sum _{a=1}^{2601}{({\mathit{ϵ}}_{ij(a)}^t-{\mathit{ϵ}}_{ij(a)})}^2}}{\sqrt{\sum _{a=1}^{2601}{({\mathit{ϵ}}_{ij(a)}^t)}^2}}\end{array}$$(48) (49) $$\begin{array}{c}\hfill \text{Disp error}=\frac{||u_i^t-u_i||}{||u_i^t||}\end{array}$$(49)

Here, $u_i^t$ and ${\mathit{ϵ}}_{ij}^t$ are the target displacement and strain components, respectively, and the $i^{\prime }s$ and $j^{\prime }s$ go from 1–2 with 1 representing the x direction, and 2 representing the y direction. The counter a runs over all nodes in the mesh.

The norms calculated with these formulas were plotted as a function of iteration number in Figure 6. For iteration zero, $\mathit{u}$ was the displacement field corrupted by noise, and $\mathit{ϵ}$ was initialized to zero strain. The plots in the figure demonstrate that the errors reduce significantly with iteration number. The exact values of the errors in the field variables after the sixth iteration are displayed in Table 1. From the table, we see that the error in $u_x$ (resp. ${\mathit{ϵ}}_{xx}$) reduced from about 50% (resp. $500\%$) to about 2%.

Figure 6. Uniaxial tension test: Error in each component of the field variables as a function of iteration number.

Figure 7. Uniaxial tension test: Total error in the field variables as a function of iteration number.

The total error in strain and displacement was computed with the following formulas:(50) $$\begin{array}{cc}\hfill \text{Total strain error}& =\sqrt{\frac{||{\mathit{ϵ}}_{xx}^{err}{||}^2+||{\mathit{ϵ}}_{yy}^{err}{||}^2+2||{\mathit{ϵ}}_{xy}^{err}{||}^2}{||{\mathit{ϵ}}_{xx}^t{||}^2+||{\mathit{ϵ}}_{yy}^t{||}^2+2||{\mathit{ϵ}}_{xy}^t{||}^2}}\hfill \end{array}$$(50) (51) $$\begin{array}{cc}\hfill \text{Total disp error}& =\sqrt{\frac{||u_x^{err}{||}^2+||u_y^{err}{||}^2}{||u_x^t{||}^2+||u_y^t{||}^2}}\hfill \end{array}$$(51)

Here, $u_i^{err}=u_i-u_i^t$ and ${\mathit{ϵ}}_{ij}^{err}={\mathit{ϵ}}_{ij}-{\mathit{ϵ}}_{ij}^t$. A plot of these total errors as a function of iteration number is shown in Figure 7, which shows that the total errors decrease with iteration.

Table 1. Uniaxial tension test: Errors in the recovered field variables after the sixth iteration, and initial errors in the input field variables.

Field variable	Error in input field (%)	Error in reconstructed field (%)
$u_x$	48.8	2.09
$u_y$	0	0.072
$||\mathit{u}||$	22.6	0.759
${\mathit{ϵ}}_{xx}$	511	1.78
${\mathit{ϵ}}_{yy}$	0	0.526
${\mathit{ϵ}}_{xy}$	1883	13.5
$||\mathit{ϵ}||$	277	1.49

3.3. Biaxial tension test

For the second verification test, we consider the same body, this time under biaxial tension; c.f. Figure 2. The stress at infinity was ${\mathit{\sigma }}_{xx}={\mathit{\sigma }}_{yy}={\mathit{\sigma }}^{\infty }=6\times {10}^{-2}$. Next, Gaussian noise was added to the lateral displacement field, and reconstructions were performed using the following parameters: $\mathit{\alpha }={10}^{-5}$, $\mathit{\beta }=10$, $\mathit{\delta }={10}^{-8}$, $n=0.5$, $T_{xx}=50$, $T_{yy}=5\times {10}^4$. These are approximately the same parameters as the last experiment, except with $T_{xx}$ significantly increased to accommodate the applied horizontal stress.

Figure 8. Biaxial tension test results. Upper row: Target noiseless strain fields from uniaxial tension test. Lower row: Strain fields recovered from SPREME. Relatively large weights on noisy lateral displacement increase the error in the recovered strains relative to those in Figure 4.

The target and reconstructed strains at the sixth iteration are shown in Figure 8. The features in the lateral and shear strains are recovered, but we see significantly more error in this test case than in Figure 4. That error is quantified in Table 2. The table shows clearly that the errors in the lateral displacement, and lateral and shear strains reduced significantly from the errors in the input field. This means that the method is able to work also for this type of deformation. The method did not work as well as it did for the case when the domain was under uniaxial tension, however. This is most likely due to the need to use relatively large weights on the noisy lateral displacements to capture the effect of significant lateral normal stress applied at infinity. Thus, the effect of the noise in that displacement component was amplified. Nevertheless, the error in ${\mathit{ϵ}}_{xx}$ reduced from about $500\%$ to about $15\%$, a very significant reduction.

Table 2. Biaxial tension test: Errors in the recovered field variables after the sixth iteration, and in the input field variables for the biaxial tension experiment with noise.

Field variable	Error in input field (%)	Error in reconstructed field (%)
$u_x$	51.2	6.63
$u_y$	0	0.103
$||\mathit{u}||$	36.6	3.83
${\mathit{ϵ}}_{xx}$	495	15.4
${\mathit{ϵ}}_{yy}$	0	1.04
${\mathit{ϵ}}_{xy}$	1370	59.3
$||\mathit{ϵ}||$	422	13.9

4. Validation

In the verification section, we demonstrated that SPREME correctly reconstructs plane stress displacement fields. Here, in the validation section, we demonstrate that the SPREME implementation of the plane stress model applies to physical systems of interest. The displacement data used for the test were measured within a gelatin phantom under uniaxial compression as described in [9]. This displacement data were analysed previously by an iterative elastic modulus inversion algorithm.[9] The results from that study, and a smoothed version of the measured displacements, will be used as a reference to compare to the results produced by SPREME.

Figure 9. Axial displacement fields. First row: Measured displacement field; second row: Displacements reconstructed by SPREME; third row: Reference displacements from an iterative inversion algorithm.

Figure 10. Lateral displacement fields. First row: Measured displacements; second row: A smoothed version of the measured displacements; third row: Displacements reconstructed by SPREME; fourth row: Reference results from an iterative inversion algorithm. We note general agreement in gross features between all rows for targets 1 and 2. For targets 3 and 4, there is disagreement between the reference and the other fields, which we attribute to inaccuracies in assumed boundary conditions used to create the reference solution.

Figure 11. Axial strain fields. First row: Measured strains; second row: Strains reconstructed by SPREME; third row: Reference results from an iterative inversion algorithm. We note qualitative and quantitative agreements between the rows; we note also the presence of boundary artefacts in the reference fields that are not present in the SPREME fields.

Figure 12. Lateral strain fields. First row: Measured strains; second row: Strains reconstructed by SPREME; third row: Reference results from an iterative inversion algorithm. We note general agreement of the main features between the SPREME and reference results; boundary noise present in reference results is not apparent in SPREME results; SPREME fields show higher dynamic range than reference fields.

Figure 13. Shear strain fields. First row: Measured strains; second row: Strains reconstructed by SPREME; third row: Reference results from an iterative inversion algorithm. We note general agreement of the main features between the SPREME and reference results; boundary noise present in reference results is not apparent in SPREME results; SPREME fields show higher dynamic range than reference fields.

Figure 14. Displacement fields measured in-vivo in the breast as described in [10]. Upper row: Fibroadenoma, a benign breast tumour (FA 3) Lower row: Invasive Ductal Carcinoma, a malignant breast tumour (IDC 2).

Figure 15. Displacement and strain images for fibroadenoma 3. We note that the reconstructed SPREME lateral displacements reproduce the gross features seen in the smoothed measurements, while the reference method does not.

Figure 16. Displacement and strain images for invasive ductal carcinoma 2. In this example, we see more consensus between the gross features of the smoothed measurements, the reconstructed SPREME lateral displacements and the reference lateral displacements.

4.1. Phantom tests

The displacement data used for the test were measured within a gelatin phantom under uniaxial compression as described in [9]. Here, for completeness of presentation, we briefly describe the data collection procedure. An agar and gelatin phantom was manufactured to have the acoustic and mechanical properties of soft tissue using methods described in [45]. The phantom was a 100-mm cube containing four spherical inclusions of diameter 10 mm. The inclusion centres are coplanar in a horizontal plane separated by 30 mm centre-to-centre distance. Each of the inclusions has a different material Young’s modulus.[9]

4.1.1. Ultrasound imaging and displacement estimation

The phantom was imaged with a Siemens SONOLINE Antares (Siemens Medical Solutions USA, Inc, Malvern, PA) clinical ultrasound system, with a linear ultrasound transducer array (Siemens VFX9-4) which was pulsed at 8.89MHz.[9] A 15- cm $\times$ 15- cm compression plate, much larger than the phantom surface, was attached to the ultrasound transducer to help generate approximately uniaxial deformation of the phantom. The phantom was first imaged before any deformation was applied, and it was imaged after every 0.5% strain up until a total strain of 20% with respect to the phantom’s height.[9] During the imaging, RF data, representing the spatial distribution of the backscattered pressure field, were recorded. A modified block matching algorithm was used to estimate the displacement field from the RF data.[46] Because the models used in this paper depend on the small strain assumption, we use only the measured deformation fields corresponding to small strain; specifically, we choose $1.5\%$ overall strain.

4.1.2. Reference results

We compare our reconstructed displacement fields to two reference fields.

The first reference displacement is simply the spatially smoothed version of the measured lateral displacement field. This smoothing is achieved by performing a local spatial average with corrections at the boundaries. A smoothed displacement field reveals its main features, but, of course, obscures detail.

A second reference displacement field is computed as the result of a second inverse problem to reconstruct the linear elastic shear modulus, as described in [10]. To that end, given a current guess of the modulus field, a forward elasticity problem is solved driven by assumed boundary conditions. The modulus field is updated to minimize the difference between measured and predicted displacements. Only axial displacements are used in the minimization.

4.1.3. Phantom results and discussion

The displacement for each target was downsampled to a 63 $\times$ 54 grid, which is the same grid that was used to generate the reference results.[9] Images of the measured, smoothed, reconstructed and reference displacements are shown in Figures 9 and 10. The parameters used for SPREME processing were: $T_{xx}={10}^{-5},T_{yy}=1,\mathit{\alpha }={10}^{-3},\mathit{\beta }=10,\mathit{\delta }={10}^{-8},n=0.5$. Eleven iterations were performed for each target, though it was observed that the strains typically converged by the 5th iteration.

We note that the reconstructed axial displacements match the measured and reference axial displacements very closely. The reconstructed lateral displacements, however, were less similar to the reference displacements. This difference is most noticeable in targets 3 and 4. It is interesting to note that the general features in the SPREME displacement images for both these targets are more similar to the smoothed measurements than are the general features in the reference displacement images. This indicates that the reference lateral displacements are in error, and this is likely due to errors in the assumed boundary conditions used to compute the reference results.²

Axial strain images (${\mathit{ϵ}}_{yy}$) for each target are shown in Figure 11, while the lateral strain fields (${\mathit{ϵ}}_{xx}$) are in Figure 12. The measured (resp. reference) strains were computed by differentiating the measured (resp. reference) displacement fields using a finite difference approximation. Three observations may be made. First, except for the lateral measured strain, all strain fields show clearly the presence of a stiff inclusion of the same size and shape, and with similar strain contrasts. The high noise level in the measured lateral displacement fields, however, is magnified through differentiation to a point where it obscures all features that might be present in the data. Second, we note boundary artefacts in the reference strain fields; these are caused by the assumed boundary conditions in the reference reconstruction method. The SPREME strain fields are relatively free of such artefacts. Third, strain contrast or dynamic range is slightly (resp. significantly) higher in the SPREME axial (resp. lateral) strain reconstructions than in the reference fields. This is perhaps due to deleterious effects of regularization used in the reference method, which tends to diminish modulus, and thus reconstructed strain, contrast. Similar observations may be made when comparing shear strains, Figure 13.

We compare the measured and reconstructed strains by estimating the contrast, and the contrast-to-noise (CNR) ratio. For the CNR, we define contrast as $|{\mathit{ϵ}}_{inc}-{\mathit{ϵ}}_{bg}|$, where ${\mathit{ϵ}}_{inc}$ is the average strain over a representative region in the inclusion, and ${\mathit{ϵ}}_{bg}$ is the average strain over a representative region in the background. The noise is estimated as the standard deviation of the strain values in the background region. For target 1, the axial strain contrast in the measured strain field is estimated to be $(8\pm 3)\times {10}^{-3}$, while that in the recovered strain field is $(8\pm 0.7)\times {10}^{-3}$. The CNR is estimated to be 2.3 for the measured strain field, and 10.4 for the recovered strain field. This confirms that the strain contrast is preserved through the processing, but the CNR is significantly improved, as clearly evident in Figure 11.

5. Application to in-vivo image data

In this section, we test SPREME with in-vivo displacement data measured from patients with breast masses. The main goal of this test is to evaluate whether SPREME is sufficiently robust to provide useful results with clinical data. The in-vivo displacement data used for the tests have been processed with an iterative inversion code in a different study described in [10]. The results from that study will be used as a reference to benchmark the results produced by SPREME. These results contain displacement fields from five biopsy-proven fibroadenomas and five biopsy proven invasive ductal carcinomas, which represent the most common forms of benign and malignant breast tumours. For the present purposes of demonstration, in this section, we present results from one of each; the reconstructed displacement and strains for the other eight cases are shown in Appendix 2. For all the cases treated, we work on the same mesh as used in [10], and focus on displacement fields corresponding to roughly 1% overall strain in order to stay in the linear range. The input measured displacement fields used are shown in Figure 14.

5.1. In-vivo results and discussion

The parameters used for SPREME processing of the clinical data were: $T_{xx}={10}^{-3},T_{yy}=1,\mathit{\alpha }=5\times {10}^{-4},\mathit{\beta }=10,\mathit{\delta }={10}^{-8},n=0.5$. The reconstructed displacement and strain fields obtained after the 11th iteration are shown in Figures 15 and 16, along with the reference fields. As with the validation study, we compare to two reference fields: the smoothed measured lateral displacement fields, and the results obtained by an iterative inversion algorithm published in [10]. Generally speaking, we see excellent agreement between all three fields in the axial displacement and strain components.

In the lateral displacement and strain fields, however, there is less consistency. The SPREME lateral displacement fields tend to reproduce the gross features seen in the smoothed lateral displacement measurements; the reference results, however, sometimes agree (c.f. Figure 16), and sometimes disagree (c.f. Figure 15) with the other two. The reference displacements seem to indicate an outward motion at the left and right side of the imaged domain in all the tumour cases (even in the eight other cases shown in Appendix 2). The smoothed measurements and the predicted displacements from SPREME, however, do not indicate this type of motion in all the tumour cases. This outward motion predicted by an iterative inversion code is probably an artefact due to the assumed boundary condition of traction free sides on the left and right sides of the imaged domain, and an implied symmetry. These assumptions may not be always valid in practice.

In contrast to the lateral displacement fields, the lateral strain fields recovered from SPREME were remarkably similar to the reference strains in all cases. This is surprising because the reconstructed lateral displacements were not always similar to the reference lateral displacements. This implies that the dominant difference between the SPREME and reference lateral displacement may be a rigid body mode.

6. Summary and conclusions

In this paper, we consider the problem of reconstructing both components of a 2D vector displacement field in a heterogenous elastic medium, given a precise measurement of one of the two components and a noisy measurement of the other. This problem is motivated by applications in quasistatic ultrasound elastography, in which technological limitations allow the very precise measurement of one component of displacement, but impede the measurement of others.

Under the assumption that the material is piecewise homogenous, we proved that the momentum equations and knowledge of one component of the displacement field determine the other displacement component up to four undetermined coefficients. We then went on to derive a variational formulation that can exploit the condition of piecewise homogeneity without other a priori knowledge of modulus distribution.

We then showed verification, validation and application of a finite element implementation from this variational formulation. The verification was performed on simulated data and showed that the method described here converges quickly, and dramatically reduces error in the noisy (uncertain) displacement component. This came at the expense of a slight increase in noise in the precise components. Nevertheless, in the verification studies conducted, the overall strain error reduced from roughly O(1) (several hundred percent) to $O({10}^{-1}\text{to}{10}^{-2})$ (a few percent). The validation studies showed that this method produces realistic displacement fields from measurements in tissue mimicking phantoms. The predicted displacements tend to be more consistent with measurements than those of a reference method in the literature. Finally, the application of SPREME to data collected in-vivo demonstrates that SPREME is sufficiently robust to work in the application domain that originally motivated its development.

Further improvements in lateral displacement estimation with SPREME are possible. Since SPREME is a post-processing method, it can be used in conjunction with the other lateral displacement estimation techniques outlined in the introduction to obtain more precise estimates of the displacement field. Since these methods produce improved measurements of the lateral displacements, relatively large weights can potentially be used on the lateral displacements when processing them with SPREME. This processing can be done efficiently because SPREME converges in very few iterations.

An area that needs further work is finding a systematic method to identify good choices of the algorithmic parameters in SPREME. Regularization parameter selection is an open and active area of research in inverse problems, and is not addressed in the present study.

The SPREME approach may be generalized to other applications, as described in [42]. The idea is that one identifies the form of the momentum equation that results for a homogeneous material property distribution. This equation is then enforced almost everywhere, as indicated in Equation (3(3) $$\begin{array}{c}\hfill \mathrm{\nabla }·\mathbf{A}=0\text{a.e.}\text{in}\mathrm{\Omega }\end{array}$$(3) ). See [42] for details.

Acknowledgements

The authors will like to thank Prof. Timothy J. Hall and members of his lab at the University of Wisconsin for sharing the phantom and in-vivo displacement data.

Additional information

Funding

The support of NIH [grant number NCI-R01CA140271] and NSF [grant number 1148124], [grant number 1148111] is greatly acknowledged.

Notes

No potential conflict of interest was reported by the authors.

1 A ‘special interface’ would be one for which there exists $[c_1,c_2,c_3,c_4]=[c_1^{\ast },c_2^{\ast },c_3^{\ast },c_4^{\ast }]\ne [0,0,0,0]$ s.t. $c_{{}_1}^{\ast }(x^2-4y^2)+c_{{}_2}^{\ast }x+c_{{}_3}^{\ast }y+c_{{}_4}^{\ast }=0$ for all (x, y) on the interface.

2 We recall that the iterative inversion code performs a forward elasticity solve at every iteration. This solve requires displacement and traction boundary conditions. The lateral direction is assumed to be traction free. The measured axial displacements were fixed as boundary conditions on all four edges of the phantom.

Appendix 1

Formulation without compatibility

Given ${\mathit{u}}_m(x,y)$, $\mathit{x}\in \mathrm{\Omega }$ and ${\mathit{ϵ}}^{{}^{(k-1)}}$, find ${\mathit{ϵ}}^k$ that minimizes:(A1) $$\begin{array}{c}\hfill \mathit{\pi }[{\mathit{ϵ}}^k]={\mathit{\pi }}_{{}_O}+{\mathit{\pi }}_{{}_R}\end{array}$$(A1)

where:(A2) $$\begin{array}{c}\hfill {\mathit{\pi }}_{{}_O}=\frac{1}{2}{\int }_{\mathrm{\Omega }}\left(\mathrm{\nabla }{\mathit{u}}_m+{(\mathrm{\nabla }{\mathit{u}}_m)}^{\text{T}}-2{\mathit{ϵ}}^k\right):\mathit{T}:\left(\mathrm{\nabla }{\mathit{u}}_m+{(\mathrm{\nabla }{\mathit{u}}_m)}^{\text{T}}-2{\mathit{ϵ}}^k\right)\mathrm{d}\mathrm{\Omega }\end{array}$$(A2)

and:(A3) $$\begin{array}{c}\hfill {\mathit{\pi }}_{{}_R}=\frac{{\mathit{\alpha }}_o}{2}{\int }_{\mathrm{\Omega }}\frac{{\left(\mathrm{\nabla }·\mathbf{A}({\mathit{ϵ}}^k)\right)}^2}{[{\left(\mathrm{\nabla }·\mathbf{A}({\mathit{ϵ}}^{{}^{(k-1)}})\right)}^2+\mathit{\delta }{]}^n}\mathrm{d}\mathrm{\Omega }.\end{array}$$(A3)

In Equation (A2(A2) $$\begin{array}{c}\hfill {\mathit{\pi }}_{{}_O}=\frac{1}{2}{\int }_{\mathrm{\Omega }}\left(\mathrm{\nabla }{\mathit{u}}_m+{(\mathrm{\nabla }{\mathit{u}}_m)}^{\text{T}}-2{\mathit{ϵ}}^k\right):\mathit{T}:\left(\mathrm{\nabla }{\mathit{u}}_m+{(\mathrm{\nabla }{\mathit{u}}_m)}^{\text{T}}-2{\mathit{ϵ}}^k\right)\mathrm{d}\mathrm{\Omega }\end{array}$$(A2) ), $\mathit{T}$ is a fourth-order weighting tensor that allows more importance to be placed on the more accurate strain component. So in ultrasound measurements, where the axial strains (${\mathit{ϵ}}_{yy}$) are typically more accurate than the lateral and shear strains (${\mathit{ϵ}}_{xx},{\mathit{ϵ}}_{xy}$), the $T_{yyyy}$ weights will be larger than $T_{xxxx}$ and $T_{xyxy}$.

Equation (A1(A1) $$\begin{array}{c}\hfill \mathit{\pi }[{\mathit{ϵ}}^k]={\mathit{\pi }}_{{}_O}+{\mathit{\pi }}_{{}_R}\end{array}$$(A1) ) was minimized with respect to $\mathit{ϵ}$ to get the weak form using the same method outlined in Section 2.2. The weak form was then discretized with finite element bilinear shape functions to get a linear system of equations. A finite element within an in-house FEM code was created to solve the linear system of equations.

Figure A1. Reconstructed strains at fifth iteration.

This formulation was tested with the perfect input displacements shown in Figure 3. The following parameters were used to run the code: $T_{yyyy}=1$, $T_{xxxx}=T_{xyxy}=1\times {10}^{-9}$, $\mathit{\alpha }=1\times {10}^{-6}$, $\mathit{\delta }=1\times {10}^{-8}$, $n=0.5$. These weights were chosen in order to simulate ultrasound measurements where the lateral and shear strains are not well known. Five iterations were performed, and the reconstructions obtained at the fifth iteration are shown in Figure A1. The lateral and shear strains look very different from the target distributions shown in Figure 4.

Table A1. Errors in the recovered strains after the fifth iteration.

Field variable	Error in reconstructed field (%)
${\mathit{ϵ}}_{xx}$	125
${\mathit{ϵ}}_{yy}$	1.41
${\mathit{ϵ}}_{xy}$	695

To quantitatively determine how accurate the reconstructions were, the errors between the reconstructed strains at the fifth iteration and target strains were computed using the formula shown in Equation (48(48) $$\begin{array}{c}\hfill \text{Strain error}=\frac{||{\mathit{ϵ}}_{ij}^t-{\mathit{ϵ}}_{ij}||}{||{\mathit{ϵ}}_{ij}^t||}\approx \frac{\sqrt{\sum _{a=1}^{2601}{({\mathit{ϵ}}_{ij(a)}^t-{\mathit{ϵ}}_{ij(a)})}^2}}{\sqrt{\sum _{a=1}^{2601}{({\mathit{ϵ}}_{ij(a)}^t)}^2}}\end{array}$$(48) ). The results of this computation are displayed in Table A1. These results show that the lateral and shear strains are not accurate.

Another method used to check the output strains from the FE code is the compatibility equations which are shown below [40]:(A4) $$\begin{array}{c}\hfill {\mathit{ϵ}}_{xx,yy}+{\mathit{ϵ}}_{yy,xx}=2{\mathit{ϵ}}_{xy,xy}.\end{array}$$(A4)

The equation used to check the compatibility condition is slightly different from the equation shown in (A4(A4) $$\begin{array}{c}\hfill {\mathit{ϵ}}_{xx,yy}+{\mathit{ϵ}}_{yy,xx}=2{\mathit{ϵ}}_{xy,xy}.\end{array}$$(A4) ). The derivation of the equation used is shown below:(A5) $$\begin{array}{cc}\hfill 2{\mathit{ϵ}}_{xy}& =u_{x,y}+u_{y,x}\hfill \end{array}$$(A5) (A6) $$\begin{array}{cc}\hfill 2{\mathit{ϵ}}_{xy}-u_{x,y}& =u_{y,x}\hfill \end{array}$$(A6) (A7) $$\begin{array}{cc}\hfill {(2{\mathit{ϵ}}_{xy}-u_{x,y})}_{,y}& =u_{y,xy}\hfill \end{array}$$(A7) (A8) $$\begin{array}{cc}\hfill {(2{\mathit{ϵ}}_{xy}-u_{x,y})}_{,y}-{\mathit{ϵ}}_{yy,x}& =\mathit{\eta }.\hfill \end{array}$$(A8)

In (A8(A8) $$\begin{array}{cc}\hfill {(2{\mathit{ϵ}}_{xy}-u_{x,y})}_{,y}-{\mathit{ϵ}}_{yy,x}& =\mathit{\eta }.\hfill \end{array}$$(A8) ), $\mathit{\eta }$ is introduced as a measure of incompatibility. Setting $\mathit{\eta }=0$ in (A8(A8) $$\begin{array}{cc}\hfill {(2{\mathit{ϵ}}_{xy}-u_{x,y})}_{,y}-{\mathit{ϵ}}_{yy,x}& =\mathit{\eta }.\hfill \end{array}$$(A8) ), and taking its derivative with respect to x, we see that it is the same as Equation (A4(A4) $$\begin{array}{c}\hfill {\mathit{ϵ}}_{xx,yy}+{\mathit{ϵ}}_{yy,xx}=2{\mathit{ϵ}}_{xy,xy}.\end{array}$$(A4) ). Therefore, by knowing $u_x,{\mathit{ϵ}}_{xy}$ and ${\mathit{ϵ}}_{yy}$, we can use Equation (A8(A8) $$\begin{array}{cc}\hfill {(2{\mathit{ϵ}}_{xy}-u_{x,y})}_{,y}-{\mathit{ϵ}}_{yy,x}& =\mathit{\eta }.\hfill \end{array}$$(A8) ) to check that the strains from the FE code satisfy the compatibility equation. $\mathit{\eta }$ was evaluated pointwise within the domain for each of the experiments performed. To get a sense of how well the compatibility equation is satisfied for each experiment, we calculate the $L^2$ norm. The formula used to evaluate the $L^2$ norm is defined below:(A9) $$\begin{array}{c}\hfill ||\mathit{\eta }||=\sqrt{{\int }_{\mathrm{\Omega }}{\mathit{\eta }}^2\mathrm{d}\mathrm{\Omega }}.\end{array}$$(A9)

The $L^2$ norm of $\mathit{\eta }$ was computed for a large range of experiments where the weights were varied, and the parameters $\mathit{\alpha }$, $\mathit{\delta }$ and n were fixed to $1\times {10}^{-6}$, $1\times {10}^{-8}$ and 0.5, respectively. The results from these experiments are shown in Table A2. The results on the table demonstrate that the compatibility equation becomes increasingly violated as $T_{xxxx}$ and $T_{xyxy}$ are reduced, and the level of violation is maximum when $T_{xxxx}=T_{xyxy}={10}^{-5}$. The mesh size for the experiment is 0.02; therefore, any value of $||\mathit{\eta }||$ smaller than this mesh size can be neglected as discretization error. When $T_{xxxx}$ and $T_{xyxy}$ become smaller than ${10}^{-3}$, then the compatibility equation becomes violated for the experiment. This means that this alternate formulation cannot be used to calculate correct strains when the weights for the lateral and shear strain components are too small. When the weights for some of the strain components are small, the FE code constrains the corresponding strains with small weights to satisfy the equilibrium equations, but the strains are not forced to come from a single-valued continuous displacement field. This formulation will therefore not work on real measured data where we would wish to use small weights to exclude the strain components calculated from the imprecisely measured displacement component.

Table A2. Compatibility results using various weights.

$T_{yyyy}$	$T_{xxxx}$	$T_{xyxy}$	$||\mathit{\eta }||$
1	1	1	0.013
1	${10}^{-1}$	${10}^{-1}$	0.013
1	${10}^{-2}$	${10}^{-2}$	0.013
1	${10}^{-3}$	${10}^{-3}$	0.017
1	${10}^{-4}$	${10}^{-4}$	0.093
1	${10}^{-5}$	${10}^{-5}$	0.127
1	${10}^{-9}$	${10}^{-9}$	0.127
1	0	0	0.127

By way of comparison, the SPREME formulation yields $||\mathit{\eta }||\approx 9\times {10}^{-3}$.

Appendix 2

Clinical data reconstructions

In this section, for completeness, we show the reconstructions obtained from processing the rest of the in-vivo data from [10]. The observations made earlier apply to these results as well. 0

Figure B1. Displacement and strain images for fibroadenoma 1.

Figure B2. Displacement and strain images for fibroadenoma 2.

Figure B3. Displacement and strain images for fibroadenoma 4.

Figure B4. Displacement and strain images for fibroadenoma 5.

Figure B5. Displacement and strain images for invasive ductal carcinoma 1.

Figure B6. Displacement and strain images for invasive ductal carcinoma 3.

Figure B7. Displacement and strain images for invasive ductal carcinoma 4.

Figure B8. Displacement and strain images for invasive ductal carcinoma 5.