Gradient Based Elastic Property Reconstruction in Digital Image Correlation Elastography.

AIM: A Conjugate Gradient implementation of the Digital Image Correlation Elastography method is presented.
METHOD: The gradient is calculated using the adjoint method, requiring only two forward solutions regardless of the number of mechanical properties to reconstruct. A power-law based multi-frequency viscoelastic model is used to relate the reconstructed mechanical properties and the Digital Image Correlation surface displacement measurements. RESULT: The method is tested against harmonic surface motion fields generated through numerical simulation and measured on a silicon phantom.
CONCLUSION: Reconstruction results show the ability of the gradient reconstruction method to detect stiff internal inclusions in simulated and phantom data. Copyright:

INTRODUCTION

The mechanical properties of breast tissue are greatly dependent on the health of the tissue. In particulars, the stiffness of cancerous breast tumors increases proportionally with the severity of cancer.[1] This difference in stiffness has led to the development of elastography, a noninvasive, image-based technique to determine the stiffness distribution within tissue.[23] Most often, the imaging modalities used for elastography are ultrasound and magnetic resonance imaging.[4]

Digital Image Correlation Elastography (DICE) is an elastography technique where the mechanical properties are reconstructed through the use of three-dimensional (3D) calibrated digital images of the breast surface under actuation. A digital image correlation (DIC) procedure is applied to these images in order to obtain the surface displacement of the breast.[45] These displacement measurements then serve as an objective against which model-based displacement calculations are optimized to reconstruct the internal stiffness distribution within the breast.[6]

In previous work, the so-called Digital Image Elasto-Tomography (DIET) method has been presented, based on stochastic minimization of the displacement objective function.[7] While initial in vivo results based on the DIET method showed promising results in detecting and localizing breast tumors,[6] the genetic algorithm that formed the base of the DIET reconstruction approach was computationally expensive and impractical for a clinical imaging modality.[78] Recent work in more direct surface-based image reconstruction approaches for the elastography method has shown that with surface data of sufficient quality, model-based reconstruction methods can directly detect the internal stiffness distributions of elastic media.[9] In addition, vastly improved data acquisition capabilities for the imaging prototype developed in the authors' laboratory has provided sufficient data to look toward more direct, gradient-based image reconstruction methods.[10] In this article, we present a Conjugate Gradient (CG) based reconstruction algorithm applied to simulated data and DIC data from the surface of a tissue-mimicking silicon phantom.

METHODS

The phantoms used in this study were made from a platinum-catalyzed soft silicone gel from Factor II Inc. (A-341) mixed with silicon fluid (Factor II-V40104) to reduce its rigidity to physiological levels.

Experimental setup

The phantoms were placed in the DICE prototype, which consists of an actuator and four pairs of stereoscopic cameras (Grasshopper 2 from Point Grey®). The 3D time-harmonic motion field of the surface was acquired using Vic-Snap and Vic-3D (from Correlated Solutions®), a DIC system allowing precise measurement of the surface displacement. The actuator is a model 2025E voice coil driver (The Modal Shop®). Generally, the DIC subset size was a square of 31 × 31 pixels with a step size (distance between two consecutive evaluation points) of 7 pixels. These parameters allow the capture of around 25,000–30,000 displacement measurements covering the entire surface of the phantom, with a resolution of 1 μm. The acquisition of the displacement was made after a startup period at the start to eliminate the transient effect. The thermal effects are small enough at the low damping rate of our phantom that it can be ignored.[11]

Digital image correlation elastography data processing

The measured displacements acquired with via Vic-3D were processed with a Python script to extract the x, y, and z positions and the u, v, and w displacements for each of the measurement points. The resulting point cloud was then processed by Hoppe's surface reconstruction method.[10] This triangulated surface is then converted into a 3D linear tetrahedral finite element (FE) mesh using the CGAL.[12] The measured data points are then projected onto the surface of the mesh at the closest surface element to serve as control points for the material property reconstruction. For each point, the FFT of the displacement time-series was calculated using the FFTpack algorithm implemented in the Scipy Python library.

Digital image correlation elastography gradient reconstruction

The forward problem is defined as: Given

The material properties: The Lamé parameters, λ and μ (μ being the shear modulus), and the density (ρ) in the spatial domain Ω

The displacement boundary conditions on the boundary ΓD

The traction vector on the boundary Γ.

Find the displacement field u that satisfies the equilibrium equation for the given material model. Using a FE formulation, this equation can be written as:

[K(θ)] u = { f },      (1)

where K(θ) is the stiffness matrix defined for materials properties θ =(λ, μ. ρ), and f is the forcing vector. The equilibrium equations for a viscoelastic isotropic compressible solid under time-harmonic motion are:

Where σ is the Cauchy stress tensor for the body, ω is the actuation pulsation, u the complex time-harmonic displacement field, u the prescribed displacement field, n the surface normal and f the traction vector. The weak form of (2) is:

A (v, u; θ ) = b (v)      (3)

With: Being the inner product and σE defined as:

σ      (5)

Combining (4) and (5) we have:

To reconstruct the mechanical properties of the phantom, we use a nonlinear inversion technique based on the (CG) method.[13] This nonlinear technique involves a computational model of the time-harmonic response of materials under external excitation (the forward problem) and estimates the spatially distributed material properties by minimizing an error function. This motion error objective function, Φ(θ), is defined as:

Where u are the measured displacements, u are the corresponding displacements calculated by FE solution (1) for the current set of mechanical properties, θ, T is a linear operator used to transform the calculated displacements to the 3D-DIC measurement space and * is the complex scalar product. For this study, the mechanical properties, θ, consist of μ and μ respectively, the shear modulus (G') and the loss modulus (G''). In order to minimize this function, the CG method is used. Starting from an initial estimate of θ, the CG update of the mechanical properties at the n-th iteration is given by:

Θ      (8)

Where pn is the search direction and αn is a step length minimizing the objective function, Φ(θ). For the nonlinear CG method used here, the search direction at the n-th iteration is given by:

p      (9)

where βk is given by the Polak-Ribère formula and gk is the gradient vector. Normally, the gradient is calculated as (9):

g = J      (10)

where J is the transpose of the Jacobian matrix of uc. Each column of J is calculated by resolving:

where K is the stiffness matrix of the forward problem. For M unknown material properties, this leads to M + 1 solutions of the forward problem to obtain the gradient vector, which is computationally expensive. The adjoint method proposed by Oberai et al.[14] uses a Lagrangian formulation to obtain a set of equations that allow the calculation of the gradient vector with only two solutions of the forward problem.

The adjoint gradient formulation begins by defining the following Lagrangian L:

L (U,W;θ) = 1/2(T[u      (12)

The differential of L is:

δL = D      (13)

Where D denotes the partial derivative operator. Setting D = 0, gives:

A (δW, U;θ) = b(δW)      (14)

Which is easily verified using (3). Setting D = 0 gives:

A (W,δU;;θ) = − (T[δu],(T[u      (15)

Combining (7) with (13), (14) and (15), we now have:

δΦ = δL = D      (16)

Since δΦ = g· δθ, we can calculate g as:

g = D      (17)

Adjoint forcing in the digital image correlation elastography system

In the DICE problem, u are displacements calculated in the FE solution space and u are displacements measured in the 3D-DIC imaging space. In general, these measurements are not located at the nodes of the FE calculation, nor do they occur at a well-distributed set of points that allow for easy interpolation to the FE space. Therefore, the linear operator, T, is used to transform the calculated displacements to the 3D-DIC measurement space. To do this, each measurement point, x, is projected onto the FE surface mesh [Figure 1].

Figure 1

Projection of the point P on the surface triangle composed of node 1, 2, and 3

This projection provides the weighting coefficients for each node of the FE surface element containing the measurement point, which allows interpolation of the calculated displacements to the projected measurement point, uͨ (xi) = T(i,:)T uͨ.

The NM x NN transformation matrix, T, is generated once during the initialization phase of the algorithm, where NM is the number of measurement points and NN are the number of nodes in the FE space. Once T is calculated, all the necessary elements to calculate the adjoint forcing vector from (15), − (T[δu], [T (uͨ [θ]) − um]), are in place. For the j element of the calculated displacement field, u, the first term in this inner product, T (δ u, represents the j column of T, such that the adjoint state solution, W, can be calculated simply as:

A (δU, W;θ) = − T      (18)

By taking advantage of the self-adjoint nature of the viscoelastic bilinear operator, A.

Multi-frequency

To improve the quality of the reconstruction, we use a multi-frequency reconstruction as it introduces more information to the poorly posed inverse problem.[15] This brings a few changes to the DICE reconstruction process. First, the data are obtained directly from a multi-frequency periodic actuation signal, usually a sum of sine waves. We then compute the FFT to isolate each individual actuation frequency. Second, the multi-frequency reconstruction algorithm works by calculating the gradient for each individual frequency, as described previously, and linearly combining them to obtain a full multi-frequency gradient. Moreover, the shear modulus is now reconstructed as a power-law, in the form of:

μ = μ      (19)

Where μ0 and α are the parameters being reconstructed, rather than μ itself. The gradient formulation then becomes:

The gradient terms for the loss modulus are calculated using the same process.

Regularization

Given the poorly posed nature of the DICE inverse problem, where internal parameters are estimated based on purely surface-based measurements, regularization terms are added to the objective function, (7), to penalize elastic property distributions with sharp variations and unreasonably high or low values. This regularization also helps avoid artifacts due to model-data mismatch around the complex, friction contact region of the actuator.

Total variation regularization

To help reduce the spatial variation of the reconstructed property distributions, θ, we introduce a total variation (TV) regularization. In our case, the TV of a parameter is defined as the integral of its squared gradient over the imaging domain.[16] This term is then added to the objective function, Φ(θ). This regularization term modifies the gradient calculation shown in Eq. 17 through the addition of the corresponding TV gradient. For the gradient term corresponding to the n material property value, θ, we have:

Where αTV is the regularization weighting, δ is a term used to ensure the differentiability of the TV operator when the property gradient is close to zero and Φn is the FE shape function that supports the nth material property value.

Tikhonov regularization

Tikhonov regularization penalizes large changes in the mechanical property values from one iteration to the next, in effect limiting the size of pn in Eq.(8). In the DICE reconstruction problem, parameter estimates near the center of the phantom are least sensitive to the surface-based displacement measurements, so the Tikhonov regularization term is scaled by the square of the distance from the central axis of the phantom (independent of the axial position) to ensure that elastic properties at the center of the phantom are only adjusted in the case of strong evidence from the surface displacement data. For the gradient term corresponding to the nth material property value, θn, we have:

Where αTK is the regularization weighting, dθn is equal to θn-θ0, β is a user-defined exponent, 2 in this case, and dist is the distance from the central axis of the phantom.

Spatial filtering

To promote smooth property distributions, we apply a Gaussian spatial filtering to the parameter distribution, θ, at each iteration. The window for this filter is set as a 5 × 5 × 5 grid surrounding each material property node.

RESULTS

Choice of the total variation coefficient

In order to correctly choose the TV weighting coefficient, α, the change in the motion error, as calculated by (7), is plotted as a function of the value of the TV coefficient, as shown in Figure 2. The optimal TV weighting corresponds to the lowest possible coefficient, which still shows the effect of the regularization, corresponding to the first point along the horizontal axis of Figure 2 where the motion error increases.

Figure 2

Variation of the error in function of the total variation coefficient, for a simulation and for an experiment

As shown in Figure 2, both for simulated and phantom data, the optimal TV coefficient lies between 10−14 and 10−12. For this study, a value of 5.10−13 was chosen for αTV.

Simulation results

To validate the DICE reconstruction method, displacement fields for homogenous and heterogeneous breast-shaped geometries were calculated by FE methods. In each case, the material properties were calculated for each frequency using a power-law in the form of Equation 19, with μ0 set to 500, and α equal to 0.45 for the shear modulus and 0.32 for the loss modulus, reflecting values measured in.[10]

For the homogeneous geometry, the properties at 60 Hz were 7.2 kPa for μR and 3.3 kPa for μI. The mono-frequency reconstruction gives a mean value of 6.1 kPa for μR (min: 5.6 kPa, max: 6.9 kPa) and 2.4 kPa for μI (min 2.2 kPa, max: 3.2 kPa), representing 84% and 72% of the true value, respectively. The multi-frequency reconstruction give a mean value of 7.2 kPa for μR (min: 6.9 kPa, max: 7.6 kPa) and 3.6 kPa for μI (min 3.3 kPa, max: 3.8 kPa), representing 100.5% and 108% of the simulated value, respectively.

For the heterogeneous geometry, the shear modulus of the inclusion was set as 20 X stiffer than the surrounding material. All other properties were identical between the two materials. While the properties of the surrounding tissue were quite accurately reconstructed (103% and 104% of the simulated value), the properties of the inclusion were incorrect, although the location of the inclusion is correctly identified. Figure 3 shows the result for mono-frequency (at 60 Hz) and multi-frequency reconstructions of the heterogeneous geometry.

Figure 3

Shear modulus reconstructions for the heterogeneous geometry simulation (left), mono-frequency reconstruction (center), and multi-frequency reconstruction (right)

Experimental results

The main objective of the DICE method is to detect and localize stiff inclusions within the imaging volume. To validate the method using measurement data, two silicon phantoms (one homogeneous and one heterogeneous) were constructed. The heterogeneous phantom consisted of a half-ellipsoid with half axes of 6.5 cm, 6.5 cm, and 7.5 cm with two inclusions (1 cm and 2 cm diameter) located at 120° from each other. The inclusions are located at a depth of 4 cm and located 1 cm and 2 cm from the phantom surface for the small inclusion and large inclusion, respectively. The two inclusions were 3D printed from ABS plastic, and as they were hollow, were visible within a standard T1 magnetic resonance image (MRI). The density of the silicon phantom material was measured at roughly 900 kg/m3 and for the reconstruction; the λ modulus was set at 150 kPa, which corresponds to a Poisson's ratio of 0.48 at a shear modulus of 6.25 kPa.

The DICE reconstruction for the homogenous phantom is shown in Figure 4. As can be seen, the stiffness distribution is largely homogeneous, with a slightly stiffer region toward the center of the phantom, where the influence of the surface-based measurements is minimal.

Figure 4

Multi-frequency reconstruction of the homogeneous phantom. Shown is the Shear modulus at 60 Hz

Figures 5 and 6 show the result of the reconstruction for the heterogeneous phantom. Figure 6 shows the comparison between the T1 MRI within a slice containing both inclusions and the corresponding slice of the multi-frequency DICE reconstruction.

Figure 5

Multi-frequency reconstruction of the heterogeneous phantom. Shown is the Shear Modulus at 60 Hz

Figure 6

T1 magnetic resonance image localization of the inclusions (left), and the corresponding slice of the digital image correlation elastography reconstruction (right), showing the Shear modulus at 60 Hz. The resolution of the magnetic resonance image is 0.8 mm per pixel, the digital image correlation elastography resolution is 3 mm per pixel

Comparison with the sweep reconstruction

The material properties for the homogeneous phantom can be compared with results obtained from a brute force sweep analysis. The sweep analysis is done by calculating the value of a motion error objective function Φ, defined as:

Where ui are the measured displacements, ui are the corresponding displacements calculated by FE solution for the current set of mechanical parameters, θ, nm is the number of measurements and * is the complex scalar product. For this study, the mechanical parameters, θ, consist of G' and G“. The parameter sweep was performed over a range from 4150 Pa to 11,500 Pa for the storage modulus and 2000 Pa to 5000 Pa for the loss modulus. The sweep analysis yields, for the multifrequency parameters, log10 (θ) and α, 3.8 and 0.42 for the shear modulus, 3.77 and 0.35 for the loss modulus. The gradient reconstruction yield 3.01 and 0.35 for the shear modulus and 3.04 and 0.36 for the loss modulus. The two reconstruction methods thus give similar values.

Exact measures of viscoelastic properties in very soft gels such as those used in these experiments are hard to obtain using traditional rheometers due to resonance effects. However, the property values obtained here are in good agreement with values obtained in previous, simplified homogeneous experiments[17] as well as phantoms studies run in other laboratories using similar gels.[11]

Second phantom

To investigate the robustness of the DICE method, we performed the same experiment previously described on a second heterogeneous phantom. This phantom was smaller in size than the first, with half axes of 4.5 cm, 4.5 cm, and 6.5 cm, and contained two inclusions located at 90° to each other. The inclusions are located at mid-height, near the phantom surface. Both inclusions are roughly oval shaped, made from cut, rigid silicone gel. Figure 7 shows their location within the phantom and the corresponding slice of the DICE reconstruction.

Figure 7

(a) Location of the inclusions within the phantom. Inclusion 1 is the stiffer of the two. (b) The corresponding slice of the digital image correlation elastography reconstruction showing the location and relative stiffness of the two inclusions

The reconstruction was performed using the same algorithm and parameters described for the 1st heterogeneous phantom [Figure 8].

Figure 8

A montage of the different slices of the 2nd phantom elasticity reconstruction

While the inclusions in the reconstruction seem to appear near the bottom of the phantom, in real life, they are located at roughly mid-height. The difference is because about 1.5 cm of the phantom is not covered by the finite element method mesh, due to a lack of spatial information on this region from the 3D DIC data.

DISCUSSION

For the homogeneous phantom, the adjoint gradient method does not find a purely homogeneous material, as the shear modulus varies from around 7 kPa in the center to around 4.3 kPa on edge. However, the change is progressive with no abrupt variation, as it is the case when inclusion is present.

For the first heterogeneous phantom, we were able to successfully identify the presence and location of the stiff inclusions, as we can see from Figure 6. However, the material properties values (shear and loss modulus) of the inclusions are not accurately identified. The inclusion, being made in plastic, is considered to have an extremely high shear modulus compared to the surrounding medium. This high rigidity corresponds to extremely long mechanical wavelengths, which cannot be accurately characterized in the small size of these inclusions.

For the second heterogeneous phantom, inclusions are well located by the reconstruction algorithm. However, an artifact from the reconstruction is present on the second line of the layer. It appears as a stiff inclusion. This is believed to be due to the actuator. Tests with a ring actuator have led to the creation of a circular pattern at the same place. It is believed that it is the compression resulting from the actuation behind the actuator that creates this pattern.

Overall, the DICE method has been proven successful to locate inclusions in two different phantoms.

CONCLUSION

We were able to successfully and accurately reconstruct the material properties using simulated displacement fields for both cases, validating the code and the model used. Using real measurements, we were able to discriminate between the homogeneous and heterogeneous cases. Moreover, for the heterogeneous phantoms, we were able to adequately locate the two inclusions. More tests need to be done to determine optimal regularization levels and to assess the robustness of the DICE method.

Financial support and sponsorship

Nil.

Conflicts of interest

There are no conflicts of interest.