

ORIGINAL ARTICLE 



Year : 2020  Volume
: 45
 Issue : 3  Page : 168174 

Gradient based elastic property reconstruction in digital image correlation elastography
Alexandre Brazy, Elijah Van Houten
Department of Mechanical Engineering, University of Sherbrooke, Sherbrooke, Quebec J1L 2R1, Canada
Date of Submission  11Nov2019 
Date of Decision  10Jan2020 
Date of Acceptance  22May2020 
Date of Web Publication  13Oct2020 
Correspondence Address: Dr. Alexandre Brazy 134A rue du Conseil, Sherbrooke, Québec J1G 1H8 Canada
Source of Support: None, Conflict of Interest: None  Check 
DOI: 10.4103/jmp.JMP_99_19
Abstract   
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 powerlaw based multifrequency 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.
Keywords: Digital image correlation elastography, elastography, viscoelastic
How to cite this article: Brazy A, Houten EV. Gradient based elastic property reconstruction in digital image correlation elastography. J Med Phys 2020;45:16874 
How to cite this URL: Brazy A, Houten EV. Gradient based elastic property reconstruction in digital image correlation elastography. J Med Phys [serial online] 2020 [cited 2020 Oct 30];45:16874. Available from: https://www.jmp.org.in/text.asp?2020/45/3/168/297955 
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, imagebased technique to determine the stiffness distribution within tissue.^{[2],[3]} 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 threedimensional (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.^{[4],[5]} These displacement measurements then serve as an objective against which modelbased displacement calculations are optimized to reconstruct the internal stiffness distribution within the breast.^{[6]}
In previous work, the socalled Digital Image ElastoTomography (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.^{[7],[8]} Recent work in more direct surfacebased image reconstruction approaches for the elastography method has shown that with surface data of sufficient quality, modelbased 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, gradientbased 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 tissuemimicking silicon phantom.
Methods   
The phantoms used in this study were made from a platinumcatalyzed soft silicone gel from Factor II Inc. (A341) mixed with silicon fluid (Factor IIV40104) 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 (Grasshopper2 from Point Grey®). The 3D timeharmonic motion field of the surface was acquired using VicSnap and Vic3D (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 Vic3D 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 timeseries 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 Γ_{N}.
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:
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 timeharmonic motion are:
Where σ_{E} is the Cauchy stress tensor for the body, ω is the actuation pulsation, u the complex timeharmonic displacement field, u_{0} the prescribed displacement field, n the surface normal and f the traction vector. The weak form of (2) is:
With: Being the inner product and σ_{E} defined as:
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 timeharmonic 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_{m} are the measured displacements, u_{c}(θ) 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 3DDIC measurement space and * is the complex scalar product. For this study, the mechanical properties, θ, consist of μ_{R} and μ_{I,} 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 nth iteration is given by:
Where p_{n} 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 nth iteration is given by:
where β_{k} is given by the PolakRibère formula and g_{k} is the gradient vector. Normally, the gradient is calculated as (9):
where J^{T} is the transpose of the Jacobian matrix of u_{c}. 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:
The differential of L is:
Where D denotes the partial derivative operator. Setting D_{w}L·δW = 0, gives:
Which is easily verified using (3). Setting D_{U}L·δU = 0 gives:
Combining (7) with (13), (14) and (15), we now have:
Since δΦ= g·
δθ, we can calculate g as:
Adjoint forcing in the digital image correlation elastography system
In the DICE problem, u_{c } are displacements calculated in the FE solution space and u_{m } are displacements measured in the 3DDIC imaging space. In general, these measurements are not located at the nodes of the FE calculation, nor do they occur at a welldistributed 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 3DDIC measurement space. To do this, each measurement point, x_{i}, 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
Click here to view 
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, .
The NMxNN 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), , are in place. For the j^{th} element of the calculated displacement field, u_{cj}, thefirst term in this inner product, T (δ u_{cj}), represents the j^{th} column of T, such that the adjoint state solution, W, can be calculated simply as:
By taking advantage of the selfadjoint nature of the viscoelastic bilinear operator, A.
Multifrequency
To improve the quality of the reconstruction, we use a multifrequency 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 multifrequency periodic actuation signal, usually a sum of sine waves. We then compute the FFT to isolate each individual actuation frequency. Second, the multifrequency reconstruction algorithm works by calculating the gradient for each individual frequency, as described previously, and linearly combining them to obtain a full multifrequency gradient. Moreover, the shear modulus is now reconstructed as a powerlaw, in the form of:
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 surfacebased 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 modeldata 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^{th} material property value, θ_{n}, 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 n^{th} 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 p
_{n} in Eq.(8). In the DICE reconstruction problem, parameter estimates near the center of the phantom are least sensitive to the surfacebased 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 n^{th} material property value, θ_{n}, we have:
Where α_{TK} is the regularization weighting, dθ_{n} is equal to θ_{n}θ_{0}, β is a userdefined 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, α_{TV}, 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 thefirst 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
Click here to view 
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 breastshaped geometries were calculated by FE methods. In each case, the material properties were calculated for each frequency using a powerlaw 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 monofrequency 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 multifrequency 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 monofrequency (at 60 Hz) and multifrequency reconstructions of the heterogeneous geometry.  Figure 3: Shear modulus reconstructions for the heterogeneous geometry simulation (left), monofrequency reconstruction (center), and multifrequency reconstruction (right)
Click here to view 
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 halfellipsoid 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/m ^{3} 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 surfacebased measurements is minimal.  Figure 4: Multifrequency reconstruction of the homogeneous phantom. Shown is the Shear modulus at 60 Hz
Click here to view 
[Figure 5] and [Figure 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 multifrequency DICE reconstruction.  Figure 5: Multifrequency reconstruction of the heterogeneous phantom. Shown is the Shear Modulus at 60 Hz
Click here to view 
 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
Click here to view 
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_{m} are the measured displacements, ui_{c}(θ) 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, log_{10}(θ) 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 thefirst, 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 midheight, 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
Click here to view 
The reconstruction was performed using the same algorithm and parameters described for the 1^{st} heterogeneous phantom [Figure 8].  Figure 8: A montage of the different slices of the 2^{nd} phantom elasticity reconstruction
Click here to view 
While the inclusions in the reconstruction seem to appear near the bottom of the phantom, in real life, they are located at roughly midheight. 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 thefirst 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.
References   
1.  Samani A, Zubovits J, Plewes D. Elastic moduli of normal and pathological human breast tissues: An inversiontechniquebased investigation of 169 samples. Phys Med Biol 2007;52:156576. 
2.  Ophir J, Céspedes I, Ponnekanti H, Yazdi Y, Li X. Elastography: A quantitative method for imaging the elasticity of biological tissues. Ultrason Imaging 1991;13:11134. 
3.  Shiina T, Nightingale KR, Palmeri ML, Hall TJ, Bamber JC, Barr RG, et al. WFUMB guidelines and recommendations for clinical use of ultrasound elastography: Part 1: Basic principles and terminology. Ultrasound Med Biol 2015;41:112647. 
4.  Schreier H, Orteu JJ, Sutton MA. Image Correlation for Shape, Motion and Deformation Measurements: Basic Concepts, Theory and Applications (1 ^{st}. ed.). Springer Publishing Company, Incorporated; 2009. 
5.  Botterill T, Lotz T, Kashif A, Chase JG. Reconstructing 3D skin surface motion for the DIET breast cancer screening system. IEEE Trans Med Imaging 2014;33:110918. 
6.  Van Houten EE, Kershaw H, Lotz T, Chase JG. Localization and detection of breast cancer tumors with digital image elastotomography. Conf Proc IEEE Eng Med Biol Soc 2012;2012:26358. 
7.  Van Houten EE, Peters A, Chase JG. Phantom elasticity reconstruction with digital image elastotomography. J Mech Behav Biomed Mater 2011;4:174154. 
8.  Peters A, Chase JG, Van Houten EE. Estimating elasticity in heterogeneous phantoms using digital image elastotomography. Med Biol Eng Comput 2009;47:6776. 
9.  Yue M, Ryan F, Suman RV, Sicheng W, Sevan G. Estimating the NonHomogeneous Elastic Modulus Distribution from Surface Deformations. International Journal of Solids and Structures 2016;83: 10.1016/j.ijsolstr.2016.01.001. 
10.  Brazy A, Van Houten EE, Homogenous Shear Modulus Reconstruction at multiple frequencies using the DICE Method, unpublished observations; 2018. 
11.  Solamen LM, McGarry MD, Tan L, Weaver JB, Paulsen KD. Phantom evaluations of nonlinear inversion MR elastography. Phys Med Biol 2018;63:145021. 
12.  Hoppe H, DeRose T, Duchamp T, McDonald J, Stuetzle W. Surface reconstruction from unorganized points. In Proceedings of the 19th annual conference on Computer graphics and interactive techniques (SIGGRAPH '92). Association for Computing Machinery, New York, NY, USA. 1992. p. 718. 
13.  Tan L, McGarry MD, Van Houten EE, Ji M, Solamen L, Weaver JB, et al. Gradientbased optimization for poroelastic and viscoelastic MR elastography. IEEE Trans Med Imaging 2017;36:23650. 
14.  Oberai AA, Gokhale NH, Feijóo GR. Solution of inverse problems in elasticity imaging using the Adjoint method, Inverse Probl; 2003. 
15.  Papazoglou S, Hirsch S, Braun J, Sack I. Multifrequency inversion in magnetic resonance elastography. Phys Med Biol 2012;57:232946. 
16.  Paulsen KD, Jiang H. Enhanced frequencydomain optical image reconstruction in tissues through totalvariation minimization. Appl Opt 1996;35:344758. 
17.  Paulsen KD, Jiang H. Enhanced frequencydomain optical image reconstruction in tissues through totalvariation minimization. Appl Opt; 1996. 
[Figure 1], [Figure 2], [Figure 3], [Figure 4], [Figure 5], [Figure 6], [Figure 7], [Figure 8]
