Displacement field estimation from OCT images utilizing speckle information with applications in quantitative elastography

For the pre-processing step, we start by taking a pair of volumetric OCT scans recorded from two sequential compression states of the sample. These data, which usually span several orders of magnitude, are commonly displayed in a base-10 logarithmic scale to enhance the visibility of lower refractions. For our purpose, we take the logarithmized data and rescale it such that its brightness lies within the interval [0, 1]. Since we consider a rotationally invariant cylindrical sample here, we search for its centre and radius through lateral segmentation and circle fitting. Using this information, the sample is aligned in the x0y plane such that its centre is located at the same point in both scans. The 'height' of the cylinder thus aligns with the z-axis, and data extending the actual height of the cylinder is cut. Figure 1 depicts a rotational MIP of OCT data pre-processed in the above way. The knowledge of the origin and height of the cylinder for both compression states provides us with an initial guess of the maximal possible axial and radial displacement which can have occurred within the sample.

The next step is the detection of bubbles within the OCT data. For this, we first apply a Gaussian filter with a standard deviation of around 0.8 to 1 to the data, in order to reduce background intensity jumps and to slightly smooth the granulated appearance of bubbles. Based on the histogram statistic of the volumetric dataset, we define a brightness threshold of 0.5%–1% of the brightest pixels, which we use to binarize the volumetric data. Each bubble in the dataset thus becomes a connected group of voxels with value 1 within a background of uniform value 0. Identifying and grouping the connected components of those binary images, we obtain the different bubbles, of which we extract their centroids (by fitting them with ellipsoids) as well as their pixel-volumes. Any thereby detected bubble which is smaller in voxel-volume than a predefined threshold (between 80 and 100 voxel-volumes in the experimental data set considered below) was discarded. See figure 7 for an illustration of the outcome of the bubble detection step described above, where the procedure was applied to the OCT images depicted in figure 1.

In the third step, we aim at estimating the motion of the detected bubbles between the two recorded volumetric images.

We now compare each of the detected bubbles in the first ('start') scan with every bubble detected in the second ('end') scan. For each candidate pair of two bubbles A and B from the 'start' and 'end' scan, respectively, we compute their voxel-volumes V_A and V_B , as well as the z-coordinates of their centroids z_A and z_B . Furthermore, we calculate the distance between their two centroids (d_AB ), the distance between the cylinder axis and the 'start' centroid (${d}_{{O}_{1}A}$), the distance between the cylinder axis and the 'end' centroid (${d}_{{O}_{2}B}$), the angle between the cylinder axis and the line connecting the 'start' and 'end' centroids (α), as well as the vertical angle in the x0y plane between the line passing through the cylinder origin and the 'start' centroid, and the line passing through the cylinder origin and the 'end' centroid (φ).

For each pair of bubbles (A, B) we check if the following basic intuitive geometric criteria are satisfied:

• (a)  
  The volume of a bubble does not change significantly in the images between different compression phases, i.e.,

  $$\begin{equation*}\left\vert {V}_{A}-{V}_{B}\right\vert {< }\varepsilon .\end{equation*}$$
• (b)  
  A bubble does not move more than the displacement (applied for deforming the sample) in axial direction d_max and more than the radial expansion of the sample in lateral direction, i.e.,

  $$\begin{align*}\hfill 0{< }{d}_{AB}& {\leqslant}{d}_{\mathrm{max}},\hfill \\ \hfill {d}_{{O}_{1}A}& {< }{d}_{{O}_{2}B}.\hfill \end{align*}$$
• (c)  
  If the sample is compressed from above, then the bubbles are mainly shifted downwards (i.e., in z-direction) between the first scan and the second, i.e.,

  $$\begin{equation*}{z}_{A}{< }{z}_{B}.\end{equation*}$$
• (d)  
  In the case of rotational invariance, the sample experiences a side bulging, which implies that a bubble moves outwards from the cylinder axis in radial direction. Hence, one expects that there is either no tangential movement φ, or that it is at most very small, i.e.,

  $$\begin{equation*}\varphi {< }{\varphi }_{\mathrm{max}}.\end{equation*}$$
• (e)  
  The possible angle α between the cylinder axis and the line passing through the starting and end points of the tracked bubbles is limited by the strain applied to the sample, i.e.,

  $$\begin{equation*}{\alpha }_{\mathrm{min}}{\leqslant}\alpha {\leqslant}{\alpha }_{\mathrm{max}}.\end{equation*}$$

Two bubbles A and B are considered to match if they satisfy all of the six inequalities stated above. Hereby, the constants ɛ, d_max, α_min, α_max, and φ_max have to be suitably chosen in dependence on the strain applied in a concrete experiment. We then get an estimate of the bubble movement by computing the vector which points from the centre of mass of bubble A to the centre of mass of the matched bubble B. Two examples of the outcome of the bubble tracking step described above are shown in figure 7 (right) and figure 8. There, the procedure was applied to the OCT images depicted in figure 1, and to the corresponding complete volumetric OCT scan, respectively.

The behaviour of the above algorithm is demonstrated in section 5 below. It is both robust and efficient/fast on OCT scans with bright sharp spot reflectors. In terms of implementation, note that it can be useful not to choose one fixed ɛ in the inequality condition above, but to separate the bubbles into subsets of 'large' and 'small' bubbles, and to use a different ɛ for each one of these groups. Furthermore, both the bubble detection step and the bubble matching step can be carried out in parallel. Note that there is typically at most one bubble in the second scan which matches a bubble in the first scan, and that not necessarily every bubble can be matched with another one, in particular, since the number of bubbles in two consecutive OCT scans is not constant. However, this is not a problem, since not all bubbles need to be matched in order to improve the reconstruction quality of the subsequent displacement field estimation.

Due to theorem 3.4, one way of computing the minimizer of the functional F is by solving the variational problem (3.14), which can for example be done via a finite element (FE) method. Let {ψ₁, ..., ψ_N } be a basis of a finite dimensional subspace V_h of ${{H}^{1}\left({\Omega}\right)}^{2}$, consisting for example of piecewise linear functions on a suitable triangulation of the domain Ω. The FE approximation u_h to the solution of (3.14) then is

$$\begin{equation*}{u}_{h}=\sum _{k=1}^{N}{z}_{k}{\psi }_{k},\end{equation*}$$

where the vector z = (z_k ) is given as the solution of the matrix-vector equation Az = y, with the N × N matrix A = (A_ik ) and the N × 1 vector y = (y_i ) defined by

$$\begin{equation*}{A}_{ik}{:=}a\left({\psi }_{k},{\psi }_{i}\right)\quad \text{and}\quad {y}_{i}{:=}b\left({\psi }_{i}\right).\end{equation*}$$

Due to the ellipticity of the bilinear form a(⋅, ⋅), the matrix A is positive definite, and thus the matrix-vector system Az = y admits a stable solution. Of course, many different versions of FE methods are applicable to the variational problem, each with its own peculiarities (see for example [9, 30, 35]). Also, for the solution of the resulting matrix-vector equation a number of different solvers are possible, from simple direct solvers to preconditioned iterative solvers. In particular, due to the underlying structure of the problem, multi-level (multi-grid) ideas are preferable, which not only lend themselves to parallelization, but can also be adapted for improving the estimation of relatively large displacement fields (see section 4.3 on multi-scale approaches).

In some cases, the minimization of F via directly solving the underlying variational problem (3.14) as discussed in the previous section is either undesirable or infeasible, for example if the size of the underlying images are too large. In this case, one can apply iterative methods for the minimization of F, which we discuss in this section.

First of all, due to the representation of F given in (3.12), we can derive the continuity and differentiability results presented in the following two propositions.

Proposition 4.1. Let assumption 3.1 hold. Then for the functional F defined in (3.5) and for all $\mathbf{u},\overline{\mathbf{u}}\in {{H}^{1}\left({\Omega}\right)}^{2}$ there holds

$$\begin{equation*}\left\vert F\left(\mathbf{u}\right)-F\left(\overline{\mathbf{u}}\right)\right\vert {\leqslant}C{\Vert}\mathbf{u}-\overline{\mathbf{u}}{{\Vert}}_{{H}^{1}\left({\Omega}\right)}\left({\Vert}\mathbf{u}-\overline{\mathbf{u}}{{\Vert}}_{{H}^{1}\left({\Omega}\right)}+{\Vert}\overline{\mathbf{u}}{{\Vert}}_{{H}^{1}\left({\Omega}\right)}+1\right),\end{equation*}$$

for a constant C > 0, and thus F is continuous on ${{H}^{1}\left({\Omega}\right)}^{2}$.

Proof. It follows from the representation (3.12) together with proposition 3.3 that

$$\begin{equation*}\begin{aligned}\hfill \left\vert F\left(\mathbf{u}\right)-F\left(\overline{\mathbf{u}}\right)\right\vert & =\left\vert \frac{1}{2}a\left(\mathbf{u},\mathbf{u}\right)-b\left(\mathbf{u}\right)-\frac{1}{2}a\left(\overline{\mathbf{u}},\overline{\mathbf{u}}\right)+b\left(\overline{\mathbf{u}}\right)\right\vert \hfill \\ \hfill & =\left\vert \frac{1}{2}a\left(\mathbf{u}-\overline{\mathbf{u}},\mathbf{u}-\overline{\mathbf{u}}\right)+a\left(\mathbf{u}-\overline{\mathbf{u}},\overline{\mathbf{u}}\right)+b\left(\overline{\mathbf{u}}-\mathbf{u}\right)\right\vert \hfill \\ \hfill & {\leqslant}C\left({\Vert}\mathbf{u}-\overline{\mathbf{u}}{{\Vert}}_{{H}^{1}\left({\Omega}\right)}^{2}+{\Vert}\overline{\mathbf{u}}{{\Vert}}_{{H}^{1}\left({\Omega}\right)}{\Vert}\mathbf{u}-\overline{\mathbf{u}}{{\Vert}}_{{H}^{1}\left({\Omega}\right)}+{\Vert}\overline{\mathbf{u}}-\mathbf{u}{{\Vert}}_{{H}^{1}\left({\Omega}\right)}\right)\hfill \\ \hfill & =C{\Vert}\mathbf{u}-\overline{\mathbf{u}}{{\Vert}}_{{H}^{1}\left({\Omega}\right)}\left({\Vert}\mathbf{u}-\overline{\mathbf{u}}{{\Vert}}_{{H}^{1}\left({\Omega}\right)}+{\Vert}\overline{\mathbf{u}}{{\Vert}}_{{H}^{1}\left({\Omega}\right)}+1\right),\hfill \end{aligned}\end{equation*}$$

from which the assertion follows. □

Proposition 4.2. Let assumption 3.1 hold. Then the functional F defined in (3.5) is twice continuously Fréchet differentiable on ${{H}^{1}\left({\Omega}\right)}^{2}$ with

$$\begin{equation}{F}^{\prime }\left(\mathbf{u}\right)\left(\mathbf{h}\right)=a\left(\mathbf{u},\mathbf{h}\right)-b\left(\mathbf{h}\right)\quad \text{and}\quad {F}^{{\prime\prime}}\left(\mathbf{u}\right)\left(\mathbf{h},\mathbf{w}\right)=a\left(\mathbf{w},\mathbf{h}\right).\end{equation} \tag{ 4.1 }$$

Proof. This is a direct consequence of the representation (3.12) of F in terms of the bilinear form a(⋅, ⋅), the linear form b(⋅), and the constant c, together with the boundedness results from proposition 3.3. □

As a consequence of the above result, we get

Corollary 4.3. Let assumption 3.1 hold. Then the functional F defined in (3.5) is strictly convex on ${{H}^{1}\left({\Omega}\right)}^{2}$.

Proof. Due to propositions 4.2 and 3.3 there holds

$$\begin{equation*}{F}^{{\prime\prime}}\left(\mathbf{u}\right)\left(\mathbf{h},\mathbf{h}\right)=a\left(\mathbf{h},\mathbf{h}\right){\geqslant}{C}_{E}{\Vert}\mathbf{h}{{\Vert}}_{{H}^{1}\left({\Omega}\right)}^{2},\end{equation*}$$

and thus it now follows from [7, proposition 17.3] that F is convex. □

The above results imply that any first- or second-order iterative method can applied for the minimization of F, especially since due to the strict convexity of F (see corollary 4.3) the minimizer is unique. In particular, the first-order optimality condition

$$\begin{equation*}{F}^{\prime }\left(\mathbf{u}\right)h=0\enspace ,\quad \forall \enspace h\in {{H}^{1}\left({\Omega}\right)}^{2},\end{equation*}$$

is both necessary and sufficient in this case (compare with theorem 3.4). Since, by the Riesz representation theorem, there exists a unique element $\nabla F\left(\mathbf{u}\right)\in {{H}^{1}\left({\Omega}\right)}^{2}$ (the gradient) such that

$$\begin{equation*}{F}^{\prime }\left(\mathbf{u}\right)h=\left\langle \enspace \nabla F\left(\mathbf{u}\right),h\enspace \right\rangle \enspace ,\quad \forall \enspace h\in {{H}^{1}\left({\Omega}\right)}^{2},\end{equation*}$$

the first-order optimality condition is equivalent to

$$\begin{equation*}\nabla F\left(\mathbf{u}\right)=0.\end{equation*}$$

Applying for example the idea of fixed-point iteration to this equation yields the simple and well-known gradient descent method

$$\begin{equation*}{\mathbf{u}}_{k+1}={\mathbf{u}}_{k}-{\omega }_{k}{F}^{\prime }\left({\mathbf{u}}_{k}\right),\end{equation*}$$

where ω_k is a suitably chosen stepsize. For our particular problem, it follows from the explicit expression of F'(u)h that the gradient ${F}^{\prime }\left(\mathbf{u}\right)\in {{H}^{1}\left({\Omega}\right)}^{2}$ is given as the unique solution of the variational problem

$$\begin{equation*}{\left\langle \enspace {F}^{\prime }\left(\mathbf{u}\right),\mathbf{v}\enspace \right\rangle }_{{H}^{1}\left({\Omega}\right)}=a\left(\mathbf{u},\mathbf{v}\right)-b\left(\mathbf{v}\right)\quad \forall \enspace \mathbf{v}\in {{H}^{1}\left({\Omega}\right)}^{2},\end{equation*}$$

and can therefore be easily computed numerically, using for example again a finite-element method, which in this case can be implemented very efficiently, since the bilinear form is now just the scalar product of ${{H}^{1}\left({\Omega}\right)}^{2}$.

Of course there is a multitude of different iterative methods which can be applied to our minimization problem. The above considerations were mainly made to underline that since we showed that the problem is well-behaved in many ways, most of them are applicable and should lead to good results.

In this section, we consider the so-called multi-scale approach (see for example [26, 31–33]), for improving the reconstruction quality of the proposed methods in a numerical implementation, in particular for larger displacement fields. As the name suggests, one considers the optical flow problem on various scales, i.e., levels of discretization, by a suitable downscaling of the input images. On each of the scales, one solves the optical flow problem and then combines the results from each scale to obtain a final optical flow estimation. The idea behind this approach is that by considering the problem on multiple scales, information about the flow on those scales can be extracted more reliably than from the original input images alone. For example, due to the inherent locality of their estimation approaches, many optical flow algorithms often end up producing flow fields on fine scales only, and are not able to capture larger displacements. The multi-scale approach is a way to produce flow fields capturing larger displacements as well.

Following the ideas and notation of [31], we start by fixing a number N of scales, indexed by s and ranging from the finest (original) scale s = 0 to the coarsest scale s = N − 1. Fixing some down-sampling factor η ∈ (0, 1), the image intensity function on the scale s, denoted by I^s (x, t), is then computed recursively via the formula

$$\begin{equation*}{I}^{s}\left(\eta \mathbf{\text{x}},t\right){:=}{G}_{\sigma }{\ast}{I}^{s-1}\left(\mathbf{\text{x}},t\right),\end{equation*}$$

where G_σ is a Gaussian filter kernel with standard deviation σ and I⁰(x, t) is the image intensity function of the original image. The value of σ depends on η and is chosen as

$$\begin{equation}\sigma \left(\eta \right){:=}{\sigma }_{0}\sqrt{{\eta }^{-2}-1},\end{equation} \tag{ 4.2 }$$

where the authors of [31] suggest σ₀ := 0.6 from experience with numerical examples. The multi-scale approach then consists in solving the optical flow problem on the scale s from the images I^s and to use the obtained displacement field u^s as an initial information for the computation on the scale s − 1. Besides a method for solving the optical flow problem on any scale s, all that is required in addition for the application of the multi-scale approach is an up-scaling routine to transfer the obtained solution u^s to the finer scale s − 1. This can for example be achieved by a suitable interpolation routine.

For the iterative solution approaches presented above, incorporating the multi-scale idea is straightforward. Starting from the coarsest scale, one applies the chosen iterative method until a suitable approximation of the displacement field is obtained. This field is then transferred to a finer grid and serves as the initial guess for the iteration on that scale, and so on.

For the direct methods discussed above, the multi-scale approach combines well with a multi-grid FE approach [9]. For this, we consider a sequence ${\left({V}_{h}^{s}\right)}_{s=0}^{N-1}$ of nested FE spaces and denote by a^s (⋅, ⋅) and b^s (⋅) the bilinear and linear forms a(⋅, ⋅) and b(⋅), in which ∇I and I_t are replaced by their Gaussian-filtered counterparts ∇I^s and ${I}_{t}^{s}$, respectively. On the coarsest level s = N − 1, one then computes the flow field ${\mathbf{u}}^{N-1}\in {V}_{h}^{N-1}$ as the solution of the variational problem

$$\begin{equation*}{a}^{N-1}\left({\mathbf{u}}^{N-1},\mathbf{v}\right)={b}^{N-1}\left(\mathbf{v}\right)\quad \forall \enspace \mathbf{v}\in {V}_{h}^{N-1}.\end{equation*}$$

For all other scales 0 ⩽ s ⩽ N − 2, the flow field ${\mathbf{u}}^{s}\in {V}_{h}^{s}$ is computed recursively via

$$\begin{equation*}{\mathbf{u}}^{s}={P}_{s-1}^{s}{\mathbf{u}}^{s-1}+{\mathbf{h}}^{s},\end{equation*}$$

where ${P}_{s-1}^{s}$ denotes the orthogonal projector from ${V}_{h}^{s-1}$ to ${V}_{h}^{s}$ and ${\mathbf{h}}^{s}\in {V}_{h}^{s}$ is defined as the (unique) solution of the variational problem

$$\begin{equation*}{a}^{s}\left({\mathbf{h}}^{s},\mathbf{v}\right)={b}^{s}\left(\mathbf{v}\right)-{a}^{s}\left({P}_{s-1}^{s}{\mathbf{u}}^{s-1},\mathbf{v}\right)\quad \forall \enspace \mathbf{v}\in {V}_{h}^{s}.\end{equation*}$$

The flow field u⁰ computed on the finest scale ${V}_{h}^{0}$ is then chosen as the approximation to the solution of the displacement field estimation problem.

In this section, we apply our displacement field estimation approach to two problems with simulated data, the first being an optical flow test problem and the second mimicking an experimental quantitative elastography setup.

5.1.1. Moving squares

The first problem, which is for example used in [3, 42] for testing optical flow estimation methods, is the problem of the moving squares. Therein, one considers two squares of constant image intensity moving towards each other. This is a difficult test problem, since the gradient of the images only contains information along the borders of the squares, and it admits a number of different plausible solutions. It is well-known that the reconstruction quality increases if for example the images are overlaid with an inhomogeneous background.

In order to test our proposed approach, we thus consider the two images I₁ and I₂ depicted in figure 2, which show two moving squares with added bubbles, as they would for example result from an OCT scan. Also depicted in figure 2 is the displacement field acting between the two images, which we want to recover.

Zoom In Zoom Out Reset image size

We want to show that utilizing the bubble motion information leads to an improved displacement field reconstruction. Hence, for the following test we assume that the displacement of all the 50 bubbles in the test images is known. However, in order to model the unavoidable inaccuracies in the tracking process, we add 0.1% relative noise to the whole set of the displacement vectors, which is already much compared to the magnitude of a single displacement vector.

For obtaining the results presented below, we use the following choice of parameters: for the smoothness-regularization, we choose α = 0.8, and when utilizing the bubble information, we choose β = 4 as well as σ = 5 in the definition of $g\left(\mathbf{\text{x}},{\hat{\mathbf{\text{x}}}}^{i}\right)$ in (3.4). For using the multi-scale approach, we take 5 scale-levels, σ₀ = 0.6, and the grid-downscaling parameter η = 0.5, which leads to σ(η) ≈ 1.03923 in formula (4.2).

Figure 3 depicts the results of the displacement field estimation, for the four different combinations of not utilizing (β = 0) or utilizing (β > 0) the bubble information, and using or not using the multi-scale approach. Both the colour and the length of the displacement vectors in the figures thereby correspond to the local magnitudes of the reconstructed displacement fields. One can clearly see that utilizing the bubble information leads to reconstructions with greatly improved quality. In particular, the general tendency of optical flow to visualize only the movement of the boundary is counteracted. Both qualitative and quantitative features of the exact displacement field (compare with figure 2)) can be recovered accurately utilizing the bubble displacement information. Furthermore, one can also see from the plots that using the multi-scale approach helps to obtain a more uniform field that is also closer in magnitude to the exact solution.

Zoom In Zoom Out Reset image size

5.1.2. Simulated elastography setup

The second test problem which we now consider is designed to mimic an experimental elastography setup. We start with the rectangular sample with a circular inclusion depicted in figure 4 (left), which we assume to be fixed on the bottom and freely moving on the sides. We apply a constant downward displacement of 10% of its size (corresponding to 20 pixels) on the top, which results in the deformed sample depicted in figure 4 (middle). The corresponding displacement field, depicted in figure 4 (right), which was used to deform the sample, was computed by solving the equations of linearized elasticity (6.2), with the Lamé parameters λ and μ as in figure 11.

Zoom In Zoom Out Reset image size

In order to simulate two OCT measurements, 200 randomly distributed bubbles were added to the sample and the resulting images before and after compression were rescaled to the interval [0, 1]. The movement of the bubbles between the images was also determined by the computed displacement field. The resulting simulated 'OCT images' together with the corresponding bubble displacement vectors ${\hat{\mathbf{u}}}^{i}$ are depicted in figure 5.

Zoom In Zoom Out Reset image size

In order to reconstruct the displacement field from the simulated OCT images, we employ our estimation procedure using the multi-scale approach. We chose the same parameter as for the moving squares example above (in particular β = 4), except for the spatial regularization parameter α, for which we now use the slightly larger choice of α = 4 compared to α = 0.8 in the previous test. The reconstructed displacement field is depicted in figure 6 (left). Note that, similarly to the exact displacement field depicted in figure 4 (right), this is only a sparse representation of the actual field, which in reality is eight times denser in each direction. For a higher resolution image of the reconstructed displacement field see figure 6 (middle). As can for example be seen from figure 6 (right), the reconstructed field agrees well with the exact field, in particular in the internal part of the sample. The total relative error is 11.53%, with the relative error in the x and y components being 9.2% and 6.96%, respectively. As the numerical examples in section 7 below show, this is sufficiently accurate to obtain quantitative material parameter reconstructions.

Zoom In Zoom Out Reset image size

After testing the proposed algorithm on simulated data, we now move to experimental data. We consider the volumetric OCT data whose rotational MIP was already depicted in figure 1 in the introduction. This data results from the consequent OCT imaging of a homogeneous sample which was compressed from above by a micrometer screw gauge. The sample itself is made from silicone rubber and is rotationally symmetric.

As a first test, we apply the bubble tracking algorithm presented in section 2 to the two OCT images shown in figure 1, the results of which are depicted in figure 7. In the left two images one can see the detected bubbles (circled in black) in the two images before and after compression. The threshold value in the algorithm was set so that only the bubbles of highest brightness would be detected. The right image depicts the resulting displacement vectors after the correspondence between the bubbles in the two images has been established by the algorithm. A visual inspection shows that the motion of the bubbles was reliably captured by our proposed tracking algorithm.

Zoom In Zoom Out Reset image size

Next, we apply our proposed algorithm to the complete 3D volumetric OCT data set, and obtain the results depicted in figure 8. Many of the brightest bubbles were identified and tracked, resulting in a set of displacement vectors, which by themselves already form a physically plausible displacement field. In fact, using the rotational symmetry of the sample, and after applying a suitable interpolation, a complete displacement field can be obtained, which is already good enough for many further applications.

Zoom In Zoom Out Reset image size

We now apply our displacement field estimation algorithm to the images depicted in figure 1, using the same settings and parameters as for the simulated data considered in the previous section, but now with 6 instead of 5 layers in the multi-scale approach, since the images are of a higher resolution. When utilizing the bubble information, we use the displacement vectors of all 672 tracked bubbles depicted in figure 8.

The resulting displacement fields are depicted in figure 9, again in the four different combinations of utilizing/not utilizing the bubble information, and using/not using the multi-scale approach. Once more it is apparent that utilizing the bubble information greatly improves the quality of the reconstructed displacement fields. While basically no correct displacement fields could be obtained with the standard Horn–Schunck method (β = 0), when including the bubble information we obtain physically meaningful fields that fit to expectations also in magnitude of the displacement. The positive influence of the multi-scale approach, although not too pronounced in this case, again shows itself in the increased uniformity of the obtained field.

Zoom In Zoom Out Reset image size

Finally, we investigate the influence of the parameter α on the reconstruction. Figure 10 depicts the obtained results for the exact same setting and parameters as before, but now for two different values of α. As one may expect, a higher value of α produces a smoother and more uniform displacement field, which is useful in physical situations where one expects the internal displacement field to be smooth anyways.

Zoom In Zoom Out Reset image size

Mathematical models linking the material parameters of a sample to the internal and external forces present in an elastography experiment are commonly derived from the equations of motion [13, 48]. Assuming that the sample is linear and isotropic, and that the deformation which the sample undergoes during the experiment is comparatively small, then the model of linearized elasticity can be used. It can for example be supplied with mixed Dirichlet and Neumann boundary conditions to model applied displacement and surface traction, respectively.

Mathematically, let Ω be a non-empty, bounded, open, connected set in ${\mathbb{R}}^{N}$, for N = 1, 2, 3, with a Lipschitz continuous boundary ∂Ω, which has two subsets Γ_D and Γ_T such that $\partial {\Omega}=\overline{{{\Gamma}}_{\mathrm{D}}\cup {{\Gamma}}_{\mathrm{T}}}$, Γ_D ∩ Γ_T = ∅ and meas(Γ_D) > 0. Given body forces f, displacement data g_D, surface traction g_T and Lamé parameters λ and μ, the (forward) problem of linearized elasticity with displacement-traction boundary conditions consists in finding the displacement field u satisfying

$$\begin{equation}\begin{aligned}\hfill -\text{div}\left(\sigma \left(\mathbf{u}\right)\right)& =\mathbf{f}\enspace ,\quad \text{in}\enspace {\Omega},\hfill \\ \hfill \mathbf{u}{\vert }_{{{\Gamma}}_{\mathrm{D}}}& ={\mathbf{g}}_{\mathrm{D}},\hfill \\ \hfill \sigma \left(\mathbf{u}\right) \overrightarrow {n}{\vert }_{{{\Gamma}}_{\mathrm{T}}}& ={\mathbf{g}}_{\mathrm{T}},\hfill \end{aligned}\end{equation} \tag{ 6.2 }$$

where $\overrightarrow {n}$ is the outward unit normal vector of ∂Ω, and the stress tensor σ defining the stress-strain relation in Ω is defined by

$$\begin{equation*}\sigma \left(\mathbf{u}\right){:=}\lambda \enspace \text{div}\left(\mathbf{u}\right)I+2\mu \mathcal{E}\left(\mathbf{u}\right),\end{equation*}$$

where I is the identity matrix and $\mathcal{E}$ is the strain tensor defined by

$$\begin{equation*}\mathcal{E}\left(\mathbf{u}\right){:=}\frac{1}{2}\left(\nabla \mathbf{u}+\nabla {\mathbf{u}}^{\mathrm{T}}\right).\end{equation*}$$

One typically considers the following weak formulation of this problem [9, 12, 19]: find a function $\mathbf{u}\in V{:=}{H}_{0,{\Gamma}}^{1}{\left({\Omega}\right)}^{N}$, where ${H}_{0,{\Gamma}}^{1}\left({\Omega}\right){:=}\left\{\mathbf{u}\in {H}^{1}\left({\Omega}\right)\enspace \vert \enspace \mathbf{u}{\vert }_{{{\Gamma}}_{\mathrm{D}}}=0\right\}$, which satisfies

$$\begin{equation}{a}_{\lambda ,\mu }\left(\mathbf{u},\mathbf{v}\right)=l\left(\mathbf{v}\right)-{a}_{\lambda ,\mu }\left({\Phi},\mathbf{v}\right)\quad \forall \mathbf{v}\in V,\end{equation} \tag{ 6.3 }$$

where Φ is a homogenization function satisfying ${\Phi}{\vert }_{{{\Gamma}}_{\mathrm{D}}}={\mathbf{g}}_{\mathrm{D}}$, and the bilinear form a_λ,μ (⋅, ⋅) and the linear form l(⋅) are given by

$$\begin{equation*}\begin{aligned}\hfill {a}_{\lambda ,\mu }\left(\mathbf{u},\mathbf{v}\right)& {:=}\underset{{\Omega}}{\int }\left(\lambda \enspace \text{div}\left(\mathbf{u}\right)\text{div}\left(\mathbf{v}\right)+2\mu \enspace \mathcal{E}\left(\mathbf{u}\right):\mathcal{E}\left(\mathbf{v}\right)\right)\enspace \enspace \mathrm{d}\mathbf{x},\hfill \\ \hfill l\left(\mathbf{v}\right)& {:=}{\left\langle \enspace \mathbf{f},\mathbf{v}\enspace \right\rangle }_{{H}^{-1}\left({\Omega}\right),{H}^{1}\left({\Omega}\right)}+{\left\langle \enspace {\mathbf{g}}_{\mathrm{T}},\mathbf{v}\enspace \right\rangle }_{{H}^{-\frac{1}{2}}\left({{\Gamma}}_{\mathrm{T}}\right),{H}^{-\frac{1}{2}}\left({{\Gamma}}_{\mathrm{T}}\right)}.\hfill \end{aligned}\end{equation*}$$

If for example f ∈ H⁻¹(Ω)^N , ${\mathbf{g}}_{\mathrm{D}}\in {H}^{\frac{1}{2}}{\left({{\Gamma}}_{\mathrm{D}}\right)}^{N}$, ${\mathbf{g}}_{\mathrm{T}}\in {H}^{-\frac{1}{2}}{\left({{\Gamma}}_{\mathrm{T}}\right)}^{N}$, and ${\Phi}\in {{H}^{1}\left({\Omega}\right)}^{N}$, then the variational problem (6.3) admits a unique solution [9, 12, 19].

After considering the forward problem of linearized elasticity in the previous section, we now turn our attention to the inverse problem, which consists in determining the Lamé parameters λ and μ from measurements of the displacement field u. Introducing for s > N/2 + 1 the parameter-to-solution map F, i.e., the nonlinear operator

$$\begin{equation*}F:\mathcal{D}\left(F\right){:=}\left\{\left(\lambda ,\mu \right)\in {H}^{s}{\left({\Omega}\right)}^{2}\enspace \vert \enspace \lambda {\geqslant}0,\mu {\geqslant}\underline{\mu }{ >}0\right\}\to {{L}^{2}\left({\Omega}\right)}^{N},\quad \left(\lambda ,\mu \right){\mapsto}\mathbf{u}\left(\lambda ,\mu \right),\end{equation*}$$

where u is the solution of the variational problem (6.3), the inverse problem can be formulated as the problem of solving the nonlinear operator equation

$$\begin{equation}F\left(\lambda ,\mu \right)=\mathbf{u}.\end{equation} \tag{ 6.4 }$$

In practical applications, instead of the true displacement field u one only has access to noisy data/measurements u^δ , which for example satisfy the estimate

$$\begin{equation}{\Vert}\mathbf{u}-{\mathbf{u}}^{\delta }{{\Vert}}_{{L}^{2}\left({\Omega}\right)}{\leqslant}\delta ,\end{equation} \tag{ 6.5 }$$

where δ ⩾ 0 denotes the noise level. Since the inverse problem is ill-posed and thus very sensitive even to a small amount of noise in the data, methods for obtaining (approximate) solutions of (6.4) need to be carefully chosen. One such choice, an iterative regularization approach based on Landweber iteration, which was used for obtaining the numerical results presented in section 7, is discussed below.

One of the most well known iterative regularization methods for solving nonlinear inverse problems is Landweber iteration [15, 21], given by

$$\begin{equation}\left({\lambda }_{k+1},{\mu }_{k+1}\right)=\left({\lambda }_{k},{\mu }_{k}\right)-{\omega }_{k}^{\delta }{F}^{\prime }{\left({\lambda }_{k},{\mu }_{k}\right)}^{{\ast}}\left(F\left({\lambda }_{k},{\mu }_{k}\right)-{\mathbf{u}}^{\delta }\right),\end{equation} \tag{ 6.6 }$$

where k is the iteration index and ${\omega }_{k}^{\delta }$ is a sequence of stepsizes. As a stopping rule one typically uses the discrepancy principle, which determines the stopping index k_* as the smallest index such that

$$\begin{equation}{\Vert}F\left({\lambda }_{k},{\mu }_{k}\right)-{\mathbf{u}}^{\delta }{\Vert}{\leqslant}\tau \delta ,\end{equation} \tag{ 6.7 }$$

where δ is as in (6.5) and τ > 1 is a constant. For the stepsizes ${\omega }_{k}^{\delta }$ one can for example use a suitably chosen constant, the minimal error, or the steepest descent stepsize [38], the latter one being given by

$$\begin{equation}{\omega }_{k}^{\delta }{:=}\frac{{\Vert}{s}_{k}^{\delta }{{\Vert}}^{2}}{{\Vert}{F}^{\prime }{\left({\lambda }_{k},{\mu }_{k}\right)}^{{\ast}}{s}_{k}^{\delta }{{\Vert}}^{2}}\enspace ,\quad {s}_{k}^{\delta }{:=}{F}^{\prime }{\left({\lambda }_{k},{\mu }_{k}\right)}^{{\ast}}\left(F\left({\lambda }_{k},{\mu }_{k}\right)-{\mathbf{u}}^{\delta }\right).\end{equation} \tag{ 6.8 }$$

In order to guarantee convergence of this, and in fact most other iterative regularization methods, one typically has to use the so-called tangential cone condition or one of its weaker variants, which are often very difficult to show for non-academic examples. Fortunately, for our present problem, if one knows the Lamé parameters in an (arbitrarily) small neighbourhood of the boundary and instead of (6.4) considers

$$\begin{equation*}{F}_{c}\left(\lambda ,\mu \right)=\mathbf{u},\end{equation*}$$

where the operator F_c incorporates the knowledge of the Lamé parameters near the boundary, then even the strong tangential cone condition holds. This guarantees the convergence of (6.6) with F replaced by F_c when combined with the discrepancy principle (6.7). For more details on this matter we refer the reader to [19].

In order to implement the iteration (6.6) and the stepsize (6.8), one needs explicit expressions for the Fréchet derivative F'(λ, μ)(h_λ , h_μ ) and its adjoint F'(λ, μ)*w. These have been given in [19] and involve the solution of a number of large-scale variational problems, which can be very costly and time consuming. Even though parallelisation can be used to reduce the computational time required for solving these variational problems, the convergence of (6.6) may still be too slow. Hence, in order to speed up the iteration, we employ Nesterov's acceleration principle, which leads to the iteration

$$\begin{equation}\begin{aligned}\hfill \left({\overline{\lambda }}_{k},{\overline{\mu }}_{k}\right)& =\left({\lambda }_{k},{\mu }_{k}\right)+{\alpha }_{k}^{\delta }\left(\left({\lambda }_{k},{\mu }_{k}\right)-\left({\lambda }_{k-1},{\mu }_{k-1}\right)\right),\hfill \\ \hfill \left({\lambda }_{k+1},{\mu }_{k+1}\right)& =\left({\overline{\lambda }}_{k},{\overline{\mu }}_{k}\right)-{\omega }_{k}^{\delta }{F}^{\prime }{\left({\overline{\lambda }}_{k},{\overline{\mu }}_{k}\right)}^{{\ast}}\left(F\left({\overline{\lambda }}_{k},{\overline{\mu }}_{k}\right)-{\mathbf{u}}^{\delta }\right),\hfill \end{aligned}\end{equation} \tag{ 6.9 }$$

where also ${\omega }_{k}^{\delta }$ is computed using the intermediate iterates $\left({\overline{\lambda }}_{k},{\overline{\mu }}_{k}\right)$. In the numerical results presented below, we used the popular choice of ${\alpha }_{k}^{\delta }{:=}\left(k-1\right)/\left(k+2\right)$; see for example [20]. For the solution of nonlinear inverse problems, iteration methods of the form (6.9) have been studied for various choices of the parameters ${\alpha }_{k}^{\delta }$ in [17, 18], and convergence was established essentially under the same tangential cone condition as needed for the convergence of Landweber iteration.

In case the noise level δ is unknown or estimates of it are unreliable, then instead of the discrepancy principle one can use heuristic stopping rules such as the heuristic discrepancy principle [23, 24], which determines the stopping index k_* by minimizing

$$\begin{equation*}k{\mapsto}\sqrt{k}{\Vert}F\left({\lambda }_{k},{\mu }_{k}\right){\Vert},\end{equation*}$$

and has been successfully employed to a wide range of both linear and nonlinear inverse problems. Even if the noise level δ is known, using the discrepancy principle can remain a difficult task due to the freedom which one has in choosing the parameter τ. Popular choices such as τ = 1.1 are often not feasible, for example when the residual cannot be brought to the level of τδ (in an acceptable amount of time) in practice. Furthermore, theoretically backed choices depend on parameters that are typically either unknown or overestimated [21]. Hence, for the numerical results presented below, we used the (heuristic) discrepancy principle together with a monitoring of the residual and the reconstruction quality as an informed guideline to manually determine a suitable stopping index for each test.

It was shown in [19, theorem 3.4] that the operator F'(λ,μ)* can be written as the composition of two operators E_s and A, where E_s basically corresponds to the adjoint of the embedding operator from H^s (Ω)² to ${{L}^{2}\left({\Omega}\right)}^{2}$. The numerical examples presented in [19] showed that dropping the operator E_s in iteration (6.6) or (6.9) leads to improved Lamé parameter reconstructions, in which larger jumps are better resolved. Mathematically, this can be understood as a (Hilbert scale) preconditioning of the problem; see for example [15]. Hence, for the numerical examples presented below, we always dropped the operator E_s in our iterative reconstruction approach (6.9).

In this section, we consider the simulated elastography setup from section 5.1.2. The results of our iterative reconstruction approach using the operator F, applied to both the exact and the estimated displacement field (see figure 6), are depicted in figure 11. Starting from the uniform initial guess (λ₀, μ₀) = (490, 10), the iteration was stopped after 85 and 27 iterations taking about 30 and 10 min, respectively, and the Young's modulus E was computed from the reconstructed Lamé parameters via (6.1). Comparing the results with figure 12 one can clearly see that the reconstructed parameters μ and E agree well with the exact parameters both in shape and in magnitude. Even though, as noted above, the relative error in the estimated displacement field is 11.53%, we still manage to obtain quantitative parameter reconstructions. In particular, the difficulty in recovering the Lamé parameter λ from just a single displacement field does not prohibit a quantitative reconstruction of the Young's modulus E.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next, we assume that the Lamé parameters λ and μ are known in a small neighbourhood of the boundary of the sample. Hence, we can use the operator F_c instead of the operator F in our iterative reconstruction approach, which we again apply to both the exact and the estimated displacement field. The results after 88 and 3 iterations taking about 30 and 1 min, respectively, are depicted in figure 13. As before, we obtain quantitative reconstructions of the Lamé parameter μ and the Young's modulus E, which are not deteriorated by the inability to reconstruct λ from only a single displacement field. Note that the obtained parameter reconstructions using the operator F_c , in particular the ones derived from the exact displacement field, better resolve the (circular) shape of the inclusion.

Zoom In Zoom Out Reset image size

After considering simulated data in the previous section, we now turn our attention to real data stemming from an OCT elastography experiment. In particular, we apply our iterative reconstruction approach to the homogeneous sample for which we have already estimated the displacement field in section 5 above.

The resulting reconstructions of the Lamé parameters λ and μ as well as the Young's modulus E, obtained after 10 iterations taking about 10 min and given in kPA, are depicted in figure 14. One can see that in the inner part of the sample a relatively uniform parameter reconstruction, in particular for μ and E, was obtained. The steep rise towards the edges of the samples is either an artefact due to the dropping of the (smoothing) operator E_s in the iterative reconstruction method, or caused by the fact that the estimated displacement field becomes somewhat inaccurate towards the sides of the sample due to the lack of detected bubbles. Computing the mean values over the inner part of the sample (pixels 300 to 700 in x-direction and 200 to 600 in y-direction), we get the parameter estimates μ = 342 kPA and E = 1015 kPA, with a standard deviation of 49 kPA and 143 kPA, respectively.

Zoom In Zoom Out Reset image size

We currently investigate the efficiency and accuracy of our displacement field estimation and material parameter reconstruction methods also on inhomogeneous samples.