|Year : 2020 | Volume
| Issue : 3 | Page : 168-174
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||11-Nov-2019|
|Date of Decision||10-Jan-2020|
|Date of Acceptance||22-May-2020|
|Date of Web Publication||13-Oct-2020|
Dr. Alexandre Brazy
134A rue du Conseil, Sherbrooke, Québec J1G 1H8
Source of Support: None, Conflict of Interest: None
| 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 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.
Keywords: Digital image correlation elastography, elastography, visco-elastic
|How to cite this article:|
Brazy A, Houten EV. Gradient based elastic property reconstruction in digital image correlation elastography. J Med Phys 2020;45:168-74
|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:168-74. 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. This difference in stiffness has led to the development of elastography, a noninvasive, image-based technique to determine the stiffness distribution within tissue., Most often, the imaging modalities used for elastography are ultrasound and magnetic resonance imaging.
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., 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.
In previous work, the so-called Digital Image Elasto-Tomography (DIET) method has been presented, based on stochastic minimization of the displacement objective function. While initial in vivo results based on the DIET method showed promising results in detecting and localizing breast tumors, the genetic algorithm that formed the base of the DIET reconstruction approach was computationally expensive and impractical for a clinical imaging modality., 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. 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. 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.
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 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.
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. This triangulated surface is then converted into a 3D linear tetrahedral finite element (FE) mesh using the CGAL. 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 Γ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 time-harmonic motion are:
Where σE is the Cauchy stress tensor for the body, ω is the actuation pulsation, u the complex time-harmonic displacement field, u0 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. 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 um are the measured displacements, uc(θ) 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 μ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 n-th iteration is given by:
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:
where βk is given by the Polak-Ribère formula and gk is the gradient vector. Normally, the gradient is calculated as (9):
where JT 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. 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 DwL·δW = 0, gives:
Which is easily verified using (3). Setting DUL·δ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, uc are displacements calculated in the FE solution space and um 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, xi, 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 jth element of the calculated displacement field, ucj, thefirst term in this inner product, T (δ ucj), represents the jth column of T, such that the adjoint state solution, W, can be calculated simply as:
By taking advantage of the self-adjoint nature of the viscoelastic bilinear operator, A.
To improve the quality of the reconstruction, we use a multi-frequency reconstruction as it introduces more information to the poorly posed inverse problem. 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:
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.
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. 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 nth 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 nth material property value.
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 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.
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.
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.
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)|
Click here to view
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/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 surface-based measurements is minimal.
|Figure 4: Multi-frequency 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 multi-frequency DICE reconstruction.
|Figure 5: Multi-frequency 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 uim are the measured displacements, uic(θ) 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  as well as phantoms studies run in other laboratories using similar gels.
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 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|
Click here to view
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|
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 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 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
Conflicts of interest
There are no conflicts of interest.
| References|| |
Samani A, Zubovits J, Plewes D. Elastic moduli of normal and pathological human breast tissues: An inversion-technique-based investigation of 169 samples. Phys Med Biol 2007;52:1565-76.
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:111-34.
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:1126-47.
Schreier H, Orteu JJ, Sutton MA. Image Correlation for Shape, Motion and Deformation Measurements: Basic Concepts, Theory and Applications (1st
. ed.). Springer Publishing Company, Incorporated; 2009.
Botterill T, Lotz T, Kashif A, Chase JG. Reconstructing 3-D skin surface motion for the DIET breast cancer screening system. IEEE Trans Med Imaging 2014;33:1109-18.
Van Houten EE, Kershaw H, Lotz T, Chase JG. Localization and detection of breast cancer tumors with digital image elasto-tomography. Conf Proc IEEE Eng Med Biol Soc 2012;2012:2635-8.
Van Houten EE, Peters A, Chase JG. Phantom elasticity reconstruction with digital image elasto-tomography. J Mech Behav Biomed Mater 2011;4:1741-54.
Peters A, Chase JG, Van Houten EE. Estimating elasticity in heterogeneous phantoms using digital image elasto-tomography. Med Biol Eng Comput 2009;47:67-76.
Yue M, Ryan F, Suman RV, Sicheng W, Sevan G. Estimating the Non-Homogeneous Elastic Modulus Distribution from Surface Deformations. International Journal of Solids and Structures 2016;83: 10.1016/j.ijsolstr.2016.01.001.
Brazy A, Van Houten EE, Homogenous Shear Modulus Reconstruction at multiple frequencies using the DICE Method, unpublished observations; 2018.
Solamen LM, McGarry MD, Tan L, Weaver JB, Paulsen KD. Phantom evaluations of nonlinear inversion MR elastography. Phys Med Biol 2018;63:145021.
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. 71-8.
Tan L, McGarry MD, Van Houten EE, Ji M, Solamen L, Weaver JB, et al
. Gradient-based optimization for poroelastic and viscoelastic MR elastography. IEEE Trans Med Imaging 2017;36:236-50.
Oberai AA, Gokhale NH, Feijóo GR. Solution of inverse problems in elasticity imaging using the Adjoint method, Inverse Probl; 2003.
Papazoglou S, Hirsch S, Braun J, Sack I. Multifrequency inversion in magnetic resonance elastography. Phys Med Biol 2012;57:2329-46.
Paulsen KD, Jiang H. Enhanced frequency-domain optical image reconstruction in tissues through total-variation minimization. Appl Opt 1996;35:3447-58.
Paulsen KD, Jiang H. Enhanced frequency-domain optical image reconstruction in tissues through total-variation minimization. Appl Opt; 1996.
[Figure 1], [Figure 2], [Figure 3], [Figure 4], [Figure 5], [Figure 6], [Figure 7], [Figure 8]