Analytical derivation of elasticity in breast phantoms for deformation tracking

Purpose Patient-specific biomedical modeling of the breast is of interest for medical applications such as image registration, image guided procedures and the alignment for biopsy or surgery purposes. The computation of elastic properties is essential to simulate deformations in a realistic way. This study presents an innovative analytical method to compute the elastic modulus and evaluate the elasticity of a breast using magnetic resonance (MRI) images of breast phantoms. Methods An analytical method for elasticity computation was developed and subsequently validated on a series of geometric shapes, and on four physical breast phantoms that are supported by a planar frame. This method can compute the elasticity of a shape directly from a set of MRI scans. For comparison, elasticity values were also computed numerically using two different simulation software packages. Results Application of the different methods on the geometric shapes shows that the analytically derived elongation differs from simulated elongation by less than 9% for cylindrical shapes, and up to 18% for other shapes that are also substantially vertically supported by a planar base. For the four physical breast phantoms, the analytically derived elasticity differs from numeric elasticity by 18% on average, which is in accordance with the difference in elongation estimation for the geometric shapes. The analytic method has shown to be multiple orders of magnitude faster than the numerical methods. Conclusion It can be concluded that the analytical elasticity computation method has good potential to supplement or replace numerical elasticity simulations in gravity-induced deformations, for shapes that are substantially supported by a planar base perpendicular to the gravitational field. The error is manageable, while the calculation procedure takes less than one second as opposed to multiple minutes with numerical methods. The results will be used in the MRI and Ultrasound Robotic Assisted Biopsy (MURAB) project.


Introduction
Screening and staging of breast cancer for diagnosis and subsequent treatment is based on medical images acquired on This project has received funding from the European Unions Horizon 2020 research and innovation programme under Grant Agreement No. 688188. tive deformation models, the elasticity of the model needs to be known with good accuracy, i.e., the difference between computed and actual elasticity must be small. In this study, we aim for a maximum difference in the order of 10%, or at most two times the elasticity variation among FEM-simulated elasticity values. Image registration techniques based on image intensities could be used for small deformations [24], but do not work in cases with large deformations such as the alignment from prone to supine configurations [4].
Deformation of the breast occurs due to body movements. Various physics-based numerical procedures have been presented for biomechanical modeling and soft tissue deformation. The most common computational schemes are based on linear or nonlinear biomechanical models including mass-spring methods (MSM) [2,7,20,23], the mass-tensor method [10,22], the boundary element method [13,17] and conventional finite element modeling (FEM) [3,25,26].
In an MSM system, an object is modeled by a collection of point masses linked together with massless springs.
Recent studies show the use of FEM to align data with large deformations of the breast [15,16]. In FEM, a body is subdivided into a set of finite elements (e.g., tetrahedral or hexahedra in 3D, triangles or other polygons in 2D). Displacements and positions of each element are approximated from discrete nodal values using interpolation functions: where h i is the interpolation function for the element containing x and φ i is the scalar weight associated with h i . Different choices for the element type and the interpolation functions exist, which depend on the accuracy requirements, geometry of the objects and computational complexity [19]. In general, FEM is used to solve a dynamic problem, which is expressed as partial differential equations (PDEs). These PDEs are then approximated with FEM. The FEM procedure has the advantage that it can handle complicated geometries (and boundaries) of high quality. A dataset of radiological 3D images of the breast anatomy (computed tomography (CT) or MRI) is required to generate a patientspecific FEM. An advantage of MRI is that it shows high sensitivity for detecting breast tumors [8]. The main FEM steps include: tissue classification/segmentation, tissue surface reconstruction, FEM volumetric mesh generation and tissue type assignment for the FEM mesh.
A patient-specific biomechanical model [11] was presented before to provide an initial deformation of the breast before registration between prone and supine MRI images. A zero-gravity reference state for both prone and supine configurations was estimated. The patient-specific unloaded configuration was obtained [12]. The biomechanical methods serve in most cases for the initialization of intensity-based image registration techniques, as in [9] or [18]. The sliding motion of the breast on the chest wall was observed [6], but usually a fixed muscle surface is applied during the FEM simulations [14,18,21].
This study introduces a method to analytically derive the elastic modulus of the breast from a pair of MRI scans, taking local differences in tissue density and elasticity into account. The two MRI scans differ by the direction of the gravitational field, which are opposite to each other. Contrary to FEMbased numerical simulations, it is not needed to convert the MRI scan into a volumetric mesh, so mechanical properties on voxel scale are preserved. Also, only one iteration over all voxels is necessary, which makes the method relatively fast The proposed analytical method requires the breast to be vertically supported by a rigid planar base. As the rib cage is approximately cylindrical, a human breast would need to be supported by a patient-mounted flat plate with a hole for the breast. In an MRI scanner, the breast coil could serve this purpose.
To avoid introduction of significant non-gravity-induced deformations when converting from prone to supine position, it is desirable to use a patient rotation system (PRS) that allows leaving the patient on the bed with breast coil attached, while being flipped over by 180°. Such a system has been developed previously by Whelan et al. [27], which theoretically could be used to take MRI scans of a planarsupported breast in both prone and supine position. It may also be possible to tilt certain MRI scanners such as the 0.25 T G-scan Brio (Esaote SpA, Genoa, Italy), although this is generally limited to rotation over 90°only.

Materials and methods
Four breast phantoms were constructed ( Fig. 1, right), consisting of a rigid base with three fiducials, stiff superficial tissue, soft deep tissue and 3-4 lesions.
The superficial and deep tissues and lesions were made of polyvinyl chloride (PVC) with plasticizer mixed in different ratios to obtain different stiffnesses. Contrary to gelatinbased phantoms, PVC is a durable material that can stay intact for extended periods. The superficial tissue consists of relatively stiff PVC which was shaped using a pair of molds ( Fig. 1, left) and afterward filled with soft PVC to mimic deep tissue. The lesions were cut in different sizes and shapes from a block of relatively stiff PVC, placed inside the deep tissue at random locations. A rigid frame was put on top and covered with a layer of stiff PVC. The four phantoms which were manufactured this way differ only in the stiffness of deep tissue and the placement of lesions. Figure 2 shows the outline of a breast phantom in a neutral reference state. Depending on the orientation (prone or supine), it is deformed by the gravitational field and tip is  displaced toward the anterior or posterior direction. The magnitude of these deformations is related to the elasticity, and the approach of the research is to reconstruct the elasticity from these deformations using different methods.
The base represents a rigid inertial frame, which must be planar and normal to the gravitational direction. While a patient's rib cage provides a rigid supportive base, it is not planar but approximately cylindrical. An external structure such as a breast coil (Fig. 2) may be required to provide this planar support.
The scanner was previously calibrated using a custom 3D calibration grid (Fig. 3, left) from which a fifth-order correction polynomial correction function was constructed. The ideal, distorted and corrected grid patterns are shown in Fig. 3. The measured residual error is 0.2 mm, so sub-pixel resolution is feasible.  The distortion-corrected MRI scans (Fig. 4, left) were segmented by intensity thresholding and automatically aligned with a rigid transformation using the three fiducials, in which the root-mean-square registration error was found to be 0.2-0.3 mm. From these data, surface and volumetric meshes in different levels of detail were constructed. Figure 4 right shows two configurations of phantom I, overlaid on each other, after segmentation and registration. A significant displacement of the tip resulting from the change in gravity field direction can be observed.

Preamble
The deformation of an object in a gravitational field is the result of elongations of tissue, which depends on the local ratio of tensile stress σ and Young's modulus E: The stress at a given location is primarily induced by the weight of the masses below that location, and also influenced by interactions with surrounding tissue. In the general case, the resulting stress distribution in the object is a complex pattern and cannot be solved analytically, requiring simulations to quantify the deformations. However, in our case, we can use the knowledge that the object's attachment to the rigid frame is planar and perpendicular to the gravity direction, when in prone and supine positions. For objects with a constant cross section such as a block or a cylinder, it can be shown (see "Analytical derivation of elasticity" We introduce the assumption that the tensile stress σ solely depends on the vertical position in the object, i.e., it is constant within any planar cross section parallel to the base. It can be shown that this assumption is valid for blocks, cylinders and prism-shaped objects which have a constant cross section. For the breast phantom shapes, the assumption can be justified by the fact that the masses of the whole breast are substantially positioned below the rigid base. To validate this assumption, the stress distribution and elongation for a range of geometric shapes are also investigated. Figure 5 schematically shows the forces and pressures acting on a shape with inhomogeneous density and elasticity, hanging from a planar, rigid attachment on the top. At a given height h, the cross-sectional area is A(h), the mass of the body below it is denoted as m(h) and the gravitational force acting on it F(h). We now derive expressions for the vertical stress σ (h) and elongation (h) for every height, leading to a formula for the displacement D of the lower extremity of the body.

Analytical derivation of elasticity
The total mass of the body up to height h is given as: The gravitational force acting on the slice at height h is calculated as: The tensile stress in the slice is generally not constant, and its exact distribution depends on many levels of tissue interactions. We are interested in the mean tensile stress σ (h), which is found by dividing the gravitational force by the slice's cross-sectional area: The tissue elasticity is also inhomogeneous in general, with local Young's modulus E(r), again averaged to E(h) for height h. The local relative elongation is = L/L 0 = σ (r)/E(r), and the mean elongation at height h is given as: The total displacement of the body's lower extremity is found by integrating all infinitesimal elongations: The purpose of this study is to find the average Young's modulus E from a pair of gravity-induced body displacements. To preserve differences in (mean) elasticity among slices, we factorize every slice's elasticity into a constant factor E and a layer-specific adjustment factorÊ(h): The displacement equation can now be written as follows: It can split into an object-specific intrinsic part which remains constant across all simulations, and an extrinsic (variable) part depending on g and E only. The intrinsic part β is defined as: Substituting into D gives: For the scanned breast phantoms, we, therefore, assume that the displacement (for small displacements) is linear in g/E, with proportionality factor β. The β value can be estimated from DICOM data, in combination with knowledge of the materials. For PVC phantoms, its density was measured to be ρ = 1.075 g/cm 3 .
Analyzing the prone and supine scans of a phantom, we have β p and β s for prone and supine, respectively. In general, β p = β s , because the shapes are significantly different: The total volume and cross-sectional area at the base are approximately equal, but due to difference in height the crosssectional shape is more squeezed in prone position than in the supine one.
The phantom height H is ill-defined due to possible irregularities at the tip, but the difference H = H p − H s can be accurately determined by comparing point clouds around the tip using, e.g., the iterative closest point algorithm [5], and optimizing H such that the total point distance is minimal, or alternatively by comparing the centroids of the point clouds.
The parameter we want to compute is the Young's modulus E. When no forces act on the phantom, it would have some shape halfway the prone and supine shapes. The tip displacement to either prone or supine shape in a gravitational field g, is H /2. We can now derive the Young's modulus E as follows:

Numerical simulation of deformations
The purpose of FEM simulations is to determine the elasticity E of the different phantoms, based on the segmented models. The general strategy is to apply a gravitational field to the FEM model of a phantom in a specific direction. This deformed model is then compared to a reference phantom which was scanned in a different orientation, providing information about the elasticity parameter. In the following subsections, we present two strategies to find the Young's modulus by simulation, of which one strategy is performed by two different simulation software packages.

Estimating theˇvalues by simulation in SOFA
In "Analytical derivation of elasticity" section, we have introduced a method to derive the values of β for the four phantoms in different orientations directly from a DICOM scan. In this section, we find β by simulation in SOFA at five different mesh resolutions [1]. For each mesh resolution, we have run a simulation with the phantom's Young's modulus set to E = 6000 Pa and gravity g = 2.0 m/s 2 . After 100 iterations, the simulation has stabilized and the vertices of the mesh in this configuration were extracted and analyzed. The displacement from the initial position follows by comparing the point clouds around the tip. The value of β then follows from Eq. (10). This procedure is repeated for each resolution of the mesh and for both prone and supine orientations, then the mean β s and β p values were computed. From the β s , β p and H , and assuming linearity of the displacement to g/E ratio, the Young's modulus E can be derived using Eqs. (11) and (12).

Supine-prone and prone-supine simulation and matching in SOFA and Febio
Taking a phantom scanned in supine configuration, the base of the phantom is immobilized and a force field sized two times the gravity (19.62 m/s 2 ) in anterior direction is applied to the phantom. After stabilization in simulation, the final state is extracted and compared to the phantom in prone position, which serves as the reference phantom.
The error value, , is defined as the distance between the simulated and reference phantoms in the area around the tip of the breast and can be positive or negative. The actual value is dependent on the elasticity parameter E of the phantom, which is optimized to bring to zero.
The minimization is performed using the Newton's method computed over E and the distance error, corrected by an adaptive step approach (when the FEM analysis software diverges). When procedure ends, i.e., when the method achieves a pre-defined error or when it reaches a maximum number of iterations, the estimated E parameter is returned with its associated error.
The procedure is then repeated for the opposite direction (prone to supine). In general, this also leads to a different E value. The mean value (square-harmonic mean-root) of E sp and E ps is then taken as the elasticity of the final phantom.

Validation of analytical stress calculation on geometric shapes
Nine homogeneous geometric shapes were generated and analyzed: two cylinders with different aspect ratios, a cone, a T-piece in normal and upside-down orientation, a half sphere, a sphere, an hourglass and a snake-like shape. Figure 6 shows the stress distribution along the vertical midway plane for all nine shapes. The first row uses the analytical computation method. The assumption that the stress distribution is constant in a cross-sectional area parallel to the base, is reflected in having constant colors in horizontal direction. The second row shows the tensile stress from numerical simulations using the SOFA software package under the same conditions. Table 1 lists the calculated and simulated β values for the same geometric shapes.
The following observations can be made:

Analytical derivation of elasticity of phantoms
Each of the four phantoms was scanned in prone and supine position, and from the resulting DICOM scans, the β p and β s values are computed using Eq. 9 and assuming a homogeneous density and elasticity distribution. From these values plus the observed vertical displacements, the E parameters are computed using Eq. 12 and the results are listed in Table 2.
It can be observed that phantom IV has the highest β and E values, making it the stiffest phantom, while phantom II is   the softest one. In general, the β values are higher in prone position, which is as expected.

Simulation ofˇin SOFA
For numerical FEM simulations, each DICOM scan was segmented and meshed at five different levels of detail and subsequently simulated in the SOFA simulation package. The resulting β values of the four phantoms (in both orientations) plus the averaged E value are listed in Table 3. Calculation of each β value requires ten simulation runs in SOFA, lasting a few minutes in total.  Figure 7 shows the analytically derived stress distribution in the transversal plane of phantom I in supine configuration together with the numerically simulated stress distribution in the same plane at low and high resolutions. It can be observed that the resulting stress patterns are comparable to that of certain geometric shapes in Fig. 6a-h. Only the analytic method shows a sharp transition at the boundary layer, as the analytical method uses slices with thickness of one voxel while the FEM-based method subdivides the volume in a different way. Table 4 lists the elasticities obtained by numerical simulation from supine to prone position and vice versa, in SOFA. As opposed to the β computation method, the prone-supine simulation method also takes nonlinearities into account which theoretically results in a more accurate estimate of the E value.

Numerical simulation by supine-prone and prone-supine matching in SOFA
For each resolution, up to ten simulation runs are needed to find the final E value in which the error vanishes. This makes the method relatively slow, requiring about twenty minutes of computation time on a quad-core 2.5 GHz computer per phantom. By parallelizing computations of the four phantoms, the total computation time for all E values was measured to be approximately half an hour.   Numerical simulation by supine-prone and prone-supine matching in FEBio  Figure 8 graphically shows the elasticity of the four phantoms, derived using the different methods, while Table 6 lists the overall phantom elasticities, averaged from the four different methods. Two methods using SOFA were presented: The first one numerically simulates the β s and β p values from the supine and prone meshes separately and measures the tip displacement D, from which the phantom's elasticity E is derived.  The second method involves finding E directly by simulation from supine to prone position such that the tip position error is eliminated. The first method seems to give consistently higher estimates for E, especially for phantom IV. Possible causes might be the nonlinearity of the displacement-tog/E ratio, i.e., β cannot be considered constant for the required range of displacements. Furthermore, the deformations of the tip resulting from proper FEM simulations influence the displacement calculations. As the second algorithm uses the iterative point cloud algorithm to minimize tip displacements and also takes nonlinear effects into account, that one can be considered more accurate than the first one.

Comparison of different elasticity measurement methods
The numerical results from FEBio simulations are in accordance with SOFA matching simulations, which is an indication that the simulations are consistent.

Discussion
We have presented a new method to analytically evaluate the elasticity of breast phantoms, from a pair of MRI scans in prone and supine position. The values found from analyzing the gravity-induced deformations are comparable to the elasticities derived from FEM simulations using FEBio and SOFA, with deviations of up to 18%. A study on nine geometric shapes has shown that the method is not only applicable to breast shapes, but also to other bodies and geometric objects as long as it is substantially supported by a planar rigid base.
The advantages of the analytical method are that the elasticity calculation is very fast (< 1 s) and takes each individually scanned voxel into account, without need for mesh generation. As the voxel intensity in a scan gives certain information about tissue type, density and/or elasticity (depending on scanning protocol), tissue inhomogeneities can be directly incorporated in the analytical computations. The main limitations are that the method is only suitable for deformations in the linear range, and that the shapes must be substantially supported by a planar base perpendicular to the gravitational field.
The fact that a human breast is relatively flexible and the chest wall is not planar but cylindrically shaped, makes clinical application difficult. An artificial planar support base could be constructed by using a patient-mounted breast coil, ideally in combination with a patient rotation system. The presented methods may also have applications in different domains, wherever deformation of bodies is involved in situations that meet the aforementioned boundary conditions The conclusion is that under specific conditions, the elasticity of a deformable object such as a human breast can be quickly computed from a pair of volumetric scans with sufficient accuracy, without need for FEM simulations. This promising result opens the door to new applications which can benefit from this complementary and near-real-time elasticity computation method.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Human and animal rights This article does not contain any studies with human participants or animals performed by any of the authors. This articles does not contain patient data.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.