Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Analytical derivation of elasticity in breast phantoms for deformation tracking

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. Keywords Biopsy · Magnetic resonance imaging · Elastic calibration · Breast Introduction different acquisition modalities and includes mammography (X-ray), ultrasound (US) and MRI. Screening and staging of breast cancer for diagnosis and sub- After image acquisition, proper localization of the tumor sequent treatment is based on medical images acquired on is essential for biopsy procedures to take tissue samples or to remove the tumor during surgery. To take full benefit from This project has received funding from the European Unions Horizon the previously acquired medical images, the location of the 2020 research and innovation programme under Grant Agreement No. tumor should be aligned from the preoperative imaging into the operating room. The position of the patient can vary from prone during MRI scanning to supine position required B Vincent Groenhuis v.groenhuis@utwente.nl for breast surgery for example. During ultrasound scanning and ultrasound-guided biopsy, the patient is returned on her Francesco Visentin francesco.visentin@univr.it back and additional compression is induced by the ultrasound probe. The computation of the elastic properties will serve University of Twente, Drienerlolaan 5, 7522 NB Enschede, as input for real-time adjustments of realistic deformations The Netherlands between preoperative and intra-operative images. For effec- University of Verona, Strada le Grazie 15, 37134 Verona, Italy 123 1642 International Journal of Computer Assisted Radiology and Surgery (2018) 13:1641–1650 tive deformation models, the elasticity of the model needs based image registration techniques, as in [9]or[18]. The to be known with good accuracy, i.e., the difference between sliding motion of the breast on the chest wall was observed computed and actual elasticity must be small. In this study, we [6], but usually a fixed muscle surface is applied during the aim for a maximum difference in the order of 10%, or at most FEM simulations [14,18,21]. two times the elasticity variation among FEM-simulated elas- This study introduces a method to analytically derive the ticity values. Image registration techniques based on image elastic modulus of the breast from a pair of MRI scans, taking intensities could be used for small deformations [24], but do local differences in tissue density and elasticity into account. not work in cases with large deformations such as the align- The two MRI scans differ by the direction of the gravitational ment from prone to supine configurations [4]. field, which are opposite to each other. Contrary to FEM- Deformation of the breast occurs due to body move- based numerical simulations, it is not needed to convert the ments. Various physics-based numerical procedures have MRI scan into a volumetric mesh, so mechanical properties been presented for biomechanical modeling and soft tissue on voxel scale are preserved. Also, only one iteration over all deformation. The most common computational schemes are voxels is necessary, which makes the method relatively fast based on linear or nonlinear biomechanical models including The proposed analytical method requires the breast to be mass-spring methods (MSM) [2,7,20,23], the mass-tensor vertically supported by a rigid planar base. As the rib cage method [10,22], the boundary element method [13,17] and is approximately cylindrical, a human breast would need to conventional finite element modeling (FEM) [3,25,26]. be supported by a patient-mounted flat plate with a hole for In an MSM system, an object is modeled by a collection the breast. In an MRI scanner, the breast coil could serve this of point masses linked together with massless springs. purpose. Recent studies show the use of FEM to align data with To avoid introduction of significant non-gravity-induced large deformations of the breast [15,16]. In FEM, a body is deformations when converting from prone to supine posi- subdivided into a set of finite elements (e.g., tetrahedral or tion, it is desirable to use a patient rotation system (PRS) hexahedra in 3D, triangles or other polygons in 2D). Dis- that allows leaving the patient on the bed with breast coil placements and positions of each element are approximated attached, while being flipped over by 180°. Such a system from discrete nodal values using interpolation functions: has been developed previously by Whelan et al. [27], which theoretically could be used to take MRI scans of a planar- supported breast in both prone and supine position. It may φ(x ) = h (x )φ (1) i i also be possible to tilt certain MRI scanners such as the 0.25 T G-scan Brio (Esaote SpA, Genoa, Italy), although this is where h is the interpolation function for the element con- generally limited to rotation over 90°only. taining x and φ is the scalar weight associated with h . i i Different choices for the element type and the interpolation functions exist, which depend on the accuracy requirements, Materials and methods geometry of the objects and computational complexity [19]. In general, FEM is used to solve a dynamic problem, which Four breast phantoms were constructed (Fig. 1, right), con- is expressed as partial differential equations (PDEs). These sisting of a rigid base with three fiducials, stiff superficial PDEs are then approximated with FEM. The FEM pro- tissue, soft deep tissue and 3–4 lesions. cedure has the advantage that it can handle complicated The superficial and deep tissues and lesions were made of geometries (and boundaries) of high quality. A dataset of polyvinyl chloride (PVC) with plasticizer mixed in different radiological 3D images of the breast anatomy (computed ratios to obtain different stiffnesses. Contrary to gelatin- tomography (CT) or MRI) is required to generate a patient- based phantoms, PVC is a durable material that can stay specific FEM. An advantage of MRI is that it shows high intact for extended periods. The superficial tissue consists of sensitivity for detecting breast tumors [8]. ThemainFEM relatively stiff PVC which was shaped using a pair of molds steps include: tissue classification/segmentation, tissue sur- (Fig. 1, left) and afterward filled with soft PVC to mimic deep face reconstruction, FEM volumetric mesh generation and tissue. The lesions were cut in different sizes and shapes from tissue type assignment for the FEM mesh. a block of relatively stiff PVC, placed inside the deep tissue at A patient-specific biomechanical model [11] was pre- random locations. A rigid frame was put on top and covered sented before to provide an initial deformation of the breast with a layer of stiff PVC. The four phantoms which were before registration between prone and supine MRI images. manufactured this way differ only in the stiffness of deep A zero-gravity reference state for both prone and supine tissue and the placement of lesions. configurations was estimated. The patient-specific unloaded Figure 2 shows the outline of a breast phantom in a neu- configuration was obtained [12]. The biomechanical meth- tral reference state. Depending on the orientation (prone or ods serve in most cases for the initialization of intensity- supine), it is deformed by the gravitational field and tip is 123 International Journal of Computer Assisted Radiology and Surgery (2018) 13:1641–1650 1643 Fig. 1 Left: pair of molds (yellow, green) for manufacturing superficial tissue (red). Right: one PVC breast phantom mounted in prone position 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. Each of the four phantoms was scanned in a 0.25 T MRI scanner (G-Scan Brio) using the 3D balanced steady- state free precession (bSSFP) sequence, with parameters TR=10ms, TE=5ms,FA = 60 , acquisition resolution 1.5 × 1.8 × 2.0 mm and isotropic reconstruction resolution 0.94 mm. Fig. 2 Breast in coil, with gravity-induced deformations in prone and supine positions (dashed lines) 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. displaced toward the anterior or posterior direction. The mag- The ideal, distorted and corrected grid patterns are shown in nitude of these deformations is related to the elasticity, and Fig. 3. The measured residual error is 0.2 mm, so sub-pixel the approach of the research is to reconstruct the elasticity resolution is feasible. from these deformations using different methods. Fig. 3 Left: MRI calibration grid. Right: actual (yellow), observed (blue) and distortion-corrected (red) grid locations of the calibration cube 123 1644 International Journal of Computer Assisted Radiology and Surgery (2018) 13:1641–1650 Fig. 4 Left: Example sagittal MRI slice. Right: Phantom I in prone and supine configuration, superimposed The distortion-corrected MRI scans (Fig. 4, left) were seg- mented 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. Elasticity estimation Preamble The deformation of an object in a gravitational field is the Fig. 5 Schematic view of force and pressure at a given height result of elongations of tissue, which depends on the local ratio of tensile stress σ and Young’s modulus E: section) that the deformation displacement can be solved analytically. The stress at a given location is primarily induced by the We introduce the assumption that the tensile stress σ weight of the masses below that location, and also influenced solely depends on the vertical position in the object, i.e., by interactions with surrounding tissue. In the general case, it is constant within any planar cross section parallel to the resulting stress distribution in the object is a complex the base. It can be shown that this assumption is valid for pattern and cannot be solved analytically, requiring simu- blocks, cylinders and prism-shaped objects which have a lations to quantify the deformations. However, in our case, constant cross section. For the breast phantom shapes, the we can use the knowledge that the object’s attachment to assumption can be justified by the fact that the masses the rigid frame is planar and perpendicular to the gravity of the whole breast are substantially positioned below the direction, when in prone and supine positions. For objects rigid base. To validate this assumption, the stress distribu- with a constant cross section such as a block or a cylinder, tion and elongation for a range of geometric shapes are also it can be shown (see “Analytical derivation of elasticity” investigated. 123 International Journal of Computer Assisted Radiology and Surgery (2018) 13:1641–1650 1645 Analytical derivation of elasticity The displacement equation can now be written as follows: g m(h) Figure 5 schematically shows the forces and pressures acting D = dh (8) on a shape with inhomogeneous density and elasticity, hang- E ˆ 0 A(h)E (h) ing from a planar, rigid attachment on the top. At a given height h, the cross-sectional area is A(h), the mass of the It can split into an object-specific intrinsic part which body below it is denoted as m(h) and the gravitational force remains constant across all simulations, and an extrinsic (variable) part depending on g and E only. The intrinsic part acting on it F (h). We now derive expressions for the vertical stress σ(h) and elongation (h) for every height, leading to β is defined as: a formula for the displacement D of the lower extremity of m(h) the body. β = dh (9) The total mass of the body up to height h is given as: ˆ 0 A(h)E (h) h Substituting into D gives: m(h) = ρ(x , y, z)dxdydz (2) D = β (10) The gravitational force acting on the slice at height h is For the scanned breast phantoms, we, therefore, assume calculated as: that the displacement (for small displacements) is linear in g/E, with proportionality factor β.The β value can be esti- F (h) = m(h)g (3) mated from DICOM data, in combination with knowledge of the materials. For PVC phantoms, its density was measured The tensile stress in the slice is generally not constant, to be ρ = 1.075 g/cm . and its exact distribution depends on many levels of tissue Analyzing the prone and supine scans of a phantom, we interactions. We are interested in the mean tensile stress σ(h), have β and β for prone and supine, respectively. In gen- p s which is found by dividing the gravitational force by the eral, β = β , because the shapes are significantly different: p s slice’s cross-sectional area: The total volume and cross-sectional area at the base are approximately equal, but due to difference in height the cross- F (h) m(h)g sectional shape is more squeezed in prone position than in σ(h) = = (4) A(h) A(h) the supine one. The phantom height H is ill-defined due to possible irreg- The tissue elasticity is also inhomogeneous in general, ularities at the tip, but the difference H = H − H can p s with local Young’s modulus E (r), again averaged to E (h) be accurately determined by comparing point clouds around for height h. The local relative elongation is  = L/L = the tip using, e.g., the iterative closest point algorithm [5], σ(r)/E (r), and the mean elongation at height h is given as: and optimizing H such that the total point distance is min- imal, or alternatively by comparing the centroids of the point σ(h) m(h)g clouds. (h) = = (5) E (h) A(h)E (h) The parameter we want to compute is the Young’s mod- ulus E. When no forces act on the phantom, it would have The total displacement of the body’s lower extremity is some shape halfway the prone and supine shapes. The tip dis- found by integrating all infinitesimal elongations: placement to either prone or supine shape in a gravitational field g,is H /2. We can now derive the Young’s modulus H H E as follows: m(h) D = (h)dh = g dh (6) A(h)E (h) 0 0 β + β p s β = (11) The purpose of this study is to find the average Young’s β g 2(β + β )g n p s E = = (12) modulus E from a pair of gravity-induced body displace- H /2 H ments. To preserve differences in (mean) elasticity among slices, we factorize every slice’s elasticity into a constant Numerical simulation of deformations factor E and a layer-specific adjustment factor E (h): The purpose of FEM simulations is to determine the elas- E (h) = E E (h) (7) ticity E of the different phantoms, based on the segmented 123 1646 International Journal of Computer Assisted Radiology and Surgery (2018) 13:1641–1650 models. The general strategy is to apply a gravitational field value. The mean value (square-harmonic mean-root) of E sp to the FEM model of a phantom in a specific direction. This and E is then taken as the elasticity of the final phantom. ps deformed model is then compared to a reference phantom which was scanned in a different orientation, providing infor- mation about the elasticity parameter. Results In the following subsections, we present two strategies to find the Young’s modulus by simulation, of which one Validation of analytical stress calculation on strategy is performed by two different simulation software geometric shapes packages. Nine homogeneous geometric shapes were generated and Estimating the ˇ values by simulation in SOFA analyzed: two cylinders with different aspect ratios, a cone, a T-piece in normal and upside-down orientation, a half sphere, In “Analytical derivation of elasticity” section, we have a sphere, an hourglass and a snake-like shape. introduced a method to derive the values of β for the four Figure 6 shows the stress distribution along the vertical phantoms in different orientations directly from a DICOM midway plane for all nine shapes. The first row uses the ana- scan. In this section, we find β by simulation in SOFA at five lytical computation method. The assumption that the stress different mesh resolutions [1]. For each mesh resolution, we distribution is constant in a cross-sectional area parallel to the have run a simulation with the phantom’s Young’s modulus base, is reflected in having constant colors in horizontal direc- set to E = 6000 Pa and gravity g = 2.0m/s . After 100 iter- tion. The second row shows the tensile stress from numerical ations, the simulation has stabilized and the vertices of the simulations using the SOFA software package under the same mesh in this configuration were extracted and analyzed. The conditions. displacement from the initial position follows by comparing Table 1 lists the calculated and simulated β values for the the point clouds around the tip. The value of β then follows same geometric shapes. from Eq. (10). This procedure is repeated for each resolution The following observations can be made: of the mesh and for both prone and supine orientations, then the mean β and β values were computed. From the β , β s p s p – For cylinder, cubic and prism-like shapes that have a and H , and assuming linearity of the displacement to g/E constant cross-sectional area (a and b), the numeri- ratio, the Young’s modulus E can be derived using Eqs. (11) cally derived stress distribution matches the analytically and (12). derived one quite well. The β values derived by both methods are well comparable (deviation under 9%). Supine–prone and prone–supine simulation and matching – For shapes that do not have a constant cross-sectional in SOFA and Febio area, but are substantially vertically supportive (c-h), the analytically calculated and SOFA-simulated β values are Taking a phantom scanned in supine configuration, the base still comparable (deviation up to 18%) although the stress of the phantom is immobilized and a force field sized two distribution is different. times the gravity (19.62 m/s ) in anterior direction is applied – For shapes in which the lower extremity is not vertically to the phantom. After stabilization in simulation, the final supported by the base, i.e., no vertical line of maximum state is extracted and compared to the phantom in prone posi- height can be drawn that entirely lies within the model tion, which serves as the reference phantom. (i), both the analytically calculated β value and the stress The error value, , is defined as the distance between the distribution are inconsistent with simulations. 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. Analytical derivation of elasticity of phantoms The minimization is performed using the Newton’s method computed over E and the distance error, corrected Each of the four phantoms was scanned in prone and supine by an adaptive step approach (when the FEM analysis soft- position, and from the resulting DICOM scans, the β and ware diverges). When procedure ends, i.e., when the method β values are computed using Eq. 9 and assuming a homoge- achieves a pre-defined error or when it reaches a maximum neous density and elasticity distribution. From these values number of iterations, the estimated E parameter is returned plus the observed vertical displacements, the E parameters with its associated error. are computed using Eq. 12 and the results are listed in Table 2. The procedure is then repeated for the opposite direction It can be observed that phantom IV has the highest β and E (prone to supine). In general, this also leads to a different E values, making it the stiffest phantom, while phantom II is 123 International Journal of Computer Assisted Radiology and Surgery (2018) 13:1641–1650 1647 Fig. 6 Analytically derived tensile stress (top row) compared with simulated stress (bottom row) for a selection of geometric shapes Table 1 Calculated and simulated β values for the nine geometric Table 3 Properties of four phantoms, derived by numerical simulation shapes in SOFA in five different resolution scales and then averaged Geometric shape Calculated β Simulated β Phantom β β HE s p a 2375 2169 I 1007 ± 58 1134 ± 38 3.28 6403 ± 207 b 2373 2229 II 947 ± 41 1125 ± 36 4.73 4297 ± 113 c 772 724 III 1131 ± 48 1259 ± 61 3.58 6549 ± 213 d 1638 1581 IV 1170 ± 43 1383 ± 34 2.93 8548 ± 184 e 4500 4979 f 213 205 g 1932 2276 Figure 7 shows the analytically derived stress distribution h 3802 3797 in the transversal plane of phantom I in supine configuration i 4942 26,499 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 cer- Table 2 Analytically derived properties of four phantoms, under the tain geometric shapes in Fig. 6a–h. Only the analytic method assumption of constant tensile stress in each cross section shows a sharp transition at the boundary layer, as the ana- Phantom β β HE s p lytical method uses slices with thickness of one voxel while the FEM-based method subdivides the volume in a different I 1215 1298 3.28 7514 way. II 1129 1269 4.73 4972 III 1356 1444 3.58 7673 IV 1420 1471 2.93 9677 Numerical simulation by supine–prone and prone–supine matching in SOFA Table 4 lists the elasticities obtained by numerical simulation the softest one. In general, the β values are higher in prone from supine to prone position and vice versa, in SOFA. As position, which is as expected. opposed to the β computation method, the prone–supine sim- ulation method also takes nonlinearities into account which Simulation of ˇ in SOFA theoretically results in a more accurate estimate of the E value. For numerical FEM simulations, each DICOM scan was For each resolution, up to ten simulation runs are needed segmented and meshed at five different levels of detail and to find the final E value in which the error vanishes. This subsequently simulated in the SOFA simulation package. The makes the method relatively slow, requiring about twenty resulting β values of the four phantoms (in both orientations) minutes of computation time on a quad-core 2.5 GHz com- plus the averaged E value are listed in Table 3. Calculation puter per phantom. By parallelizing computations of the four of each β value requires ten simulation runs in SOFA, lasting phantoms, the total computation time for all E values was a few minutes in total. measured to be approximately half an hour. 123 1648 International Journal of Computer Assisted Radiology and Surgery (2018) 13:1641–1650 Fig. 7 Tensile stress for phantom I in the transversal plane, in supine position. Left: derived using analytical method. Center and right: numerically simulated using SOFA in low resolution (center) and high resolution (right). The dashed line indicates the boundary plane between the rigid and deformable parts Table 4 Elasticity values found by numerical simulations from supine- to-prone (E ) and prone-to-supine (E ) in four different resolution sp ps scales and then averaged, using SOFA Phantom E E Mean E sp ps I 5047 ± 374 6459 ± 373 5688 ± 272 II 3395 ± 189 4513 ± 272 3895 ± 159 III 5381 ± 376 6828 ± 438 6040 ± 288 IV 6245 ± 433 7445 ± 322 6805 ± 283 Table 5 Elasticity values found by simulating from supine-to-prone (E ) and prone-to-supine (E ) in four different resolution scales and sp ps then averaged, using FEBio as software package Phantom E E Mean E Fig. 8 Young’s modulus for four phantoms, derived by four different sp ps methods I 5046 ± 272 5252 ± 307 5145 ± 254 II 4290 ± 351 4298 ± 273 4291 ± 276 III 5291.52 ± 383 5639 ± 456 5459 ± 400 Table 6 Mean elasticity values Phantom Mean E for each phantom, taken as the IV 7916 ± 1165 7564 ± 957 7731 ± 1016 average of the separate values I 6188 ± 886 derived by the four different II 4364 ± 387 methods III 6430 ± 815 IV 8190 ± 1057 Numerical simulation by supine–prone and prone–supine matching in FEBio Table 5 lists the elasticity values using the FEBio software package. The resulting elasticity values are comparable to The second method involves finding E directly by simulation those obtained by SOFA. A relatively high variance is present from supine to prone position such that the tip position error in phantom IV, which may be caused by side effects in the is eliminated. The first method seems to give consistently software package. higher estimates for E, especially for phantom IV. Possible causes might be the nonlinearity of the displacement-to- Comparison of different elasticity measurement g/E ratio, i.e., β cannot be considered constant for the methods required range of displacements. Furthermore, the defor- mations of the tip resulting from proper FEM simulations Figure 8 graphically shows the elasticity of the four phan- influence the displacement calculations. As the second algo- toms, derived using the different methods, while Table 6 lists rithm uses the iterative point cloud algorithm to minimize tip the overall phantom elasticities, averaged from the four dif- displacements and also takes nonlinear effects into account, ferent methods. that one can be considered more accurate than the first Two methods using SOFA were presented: The first one one. numerically simulates the β and β values from the supine The numerical results from FEBio simulations are in s p and prone meshes separately and measures the tip displace- accordance with SOFA matching simulations, which is an ment D, from which the phantom’s elasticity E is derived. indication that the simulations are consistent. 123 International Journal of Computer Assisted Radiology and Surgery (2018) 13:1641–1650 1649 Discussion References 1. Allard J, Cotin S, Faure F, Bensoussan PJ, Poyer F, Duriez C, We have presented a new method to analytically evaluate the Delingette H, Grisoni L (2007) Sofa-an open source framework for elasticity of breast phantoms, from a pair of MRI scans in medical simulation. In: MMVR 15-medicine meets virtual reality, prone and supine position. The values found from analyz- vol 125, pp 13–18. IOP Press ing the gravity-induced deformations are comparable to the 2. Altomonte M, Zerbato D, Botturi D, Fiorini P (2008) Simula- tion of deformable environment with haptic feedback on GPU. elasticities derived from FEM simulations using FEBio and In: IEEE/RSJ international conference on intelligent robots and SOFA, with deviations of up to 18%. A study on nine geomet- systems, 2008. IROS 2008, pp 3959–3964. IEEE ric shapes has shown that the method is not only applicable to 3. Azar FS, Metaxas DN, Schnall MD (2001) A deformable finite breast shapes, but also to other bodies and geometric objects element model of the breast for predicting mechanical deformations under external perturbations. Acad Radiol 8(10):965–975 as long as it is substantially supported by a planar rigid base. 4. Behrenbruch C, Marias K, Armitage P, Moore N, Clarke J, Brady The advantages of the analytical method are that the J (2001) Prone-supine breast MRI registration for surgical visual- elasticity calculation is very fast (< 1 s) and takes each isation. In: Medical image understanding and analysis individually scanned voxel into account, without need for 5. Besl PJ, McKay ND (1992) A method for registration of 3-D shapes. IEEE Trans Pattern Anal Mach Intell 14(2):239–256 mesh generation. As the voxel intensity in a scan gives cer- 6. Carter T, Tanner C, Beechey-Newman N, Barratt D, Hawkes D tain information about tissue type, density and/or elasticity (2008) MR navigated breast surgery: method and initial clinical (depending on scanning protocol), tissue inhomogeneities experience. In: Medical image computing and computer-assisted can be directly incorporated in the analytical computations. intervention-MICCAI 2008, pp 356–363 7. Chang YH, Chen YT, Chang CW, Lin CL (2010) Development The main limitations are that the method is only suitable for scheme of haptic-based system for interactive deformable simula- deformations in the linear range, and that the shapes must be tion. Comput Aided Des 42(5):414–424 substantially supported by a planar base perpendicular to the 8. Chevrier MC, David J, Khoury ME, Lalonde L, Labelle M, Trop gravitational field. I (2016) Breast biopsies under magnetic resonance imaging guid- ance: challenges of an essential but imperfect technique. Curr Probl The fact that a human breast is relatively flexible and the Diagn Radiol 45(3):193–204. https://doi.org/10.1067/j.cpradiol. chest wall is not planar but cylindrically shaped, makes clin- 2015.07.002 ical application difficult. An artificial planar support base 9. Conley RH, Meszoely IM, Weis JA, Pheiffer TS, Arlinghaus LR, could be constructed by using a patient-mounted breast coil, Yankeelov TE, Miga MI (2015) Realization of a biomechanical model-assisted image guidance system for breast cancer surgery ideally in combination with a patient rotation system. The using supine MRI. Int J Comput Assist Radiol Surg 10(12):1985– presented methods may also have applications in different domains, wherever deformation of bodies is involved in sit- 10. Cotin S, Delingette H, Ayache N (2000) A hybrid elastic model uations that meet the aforementioned boundary conditions for real-time cutting, deformations, and force feedback for surgery training and simulation. Vis Comput 16(8):437–452 The conclusion is that under specific conditions, the elas- 11. Eiben B, Han L, Hipwell J, Mertzanidou T, Kabus S, Bülow T, ticity of a deformable object such as a human breast can Lorenz C, Newstead G, Abe H, Keshtgar M, Ourselin S, Hawkes DJ be quickly computed from a pair of volumetric scans with (2013) Biomechanically guided prone-to-supine image registration sufficient accuracy, without need for FEM simulations. This of breast MRI using an estimated reference state. In: 2013 IEEE 10th international symposium on biomedical imaging (ISBI), pp promising result opens the door to new applications which 214–217. IEEE can benefit from this complementary and near-real-time elas- 12. Eiben B, Vavourakis V, Hipwell JH, Kabus S, Lorenz C, Buelow ticity computation method. T, Hawkes DJ (2014) Breast deformation modeling: comparison of methods to obtain a patient specific unloaded configuration. In: Proceedings of SPIE, vol 9036, pp 903615–903618 Compliance with ethical standards 13. Greminger MA, Nelson BJ (2003) Deformable object tracking using the boundary element method. In: 2003 IEEE computer soci- ety conference on computer vision and pattern recognition, 2003. Proceedings, vol 1, pp I–I. IEEE Conflict of interest The authors declare that they have no conflict of 14. Han L, Hipwell J, Mertzanidou T, Carter T, Modat M, Ourselin S, interest. Hawkes D (2011) A hybrid fem-based method for aligning prone and supine images for image guided breast surgery. In: 2011 IEEE Human and animal rights This article does not contain any studies with international symposium on biomedical imaging: from nano to human participants or animals performed by any of the authors. This macro, pp 1239–1242. IEEE articles does not contain patient data. 15. Han L, Hipwell JH, Eiben B, Barratt D, Modat M, Ourselin S, Hawkes DJ (2014) A nonlinear biomechanical model based regis- Open Access This article is distributed under the terms of the Creative tration method for aligning prone and supine MR breast images. Commons Attribution 4.0 International License (http://creativecomm IEEE Trans Med Imaging 33(3):682–694 ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, 16. Han L, Hipwell JH, Tanner C, Taylor Z, Mertzanidou T, Cardoso J, and reproduction in any medium, provided you give appropriate credit Ourselin S, Hawkes DJ (2011) Development of patient-specific to the original author(s) and the source, provide a link to the Creative biomechanical models for predicting large breast deformation. Commons license, and indicate if changes were made. Phys Med Biol 57(2):455 123 1650 International Journal of Computer Assisted Radiology and Surgery (2018) 13:1641–1650 17. James DL, Pai DK (2005) A unified treatment of elastostatic con- 24. Rueckert D, Sonoda LI, Hayes C, Hill DL, Leach MO, Hawkes tact simulation for real time haptics. In: ACM SIGGRAPH 2005 DJ (1999) Nonrigid registration using free-form deformations: courses, p 141. ACM application to breast MR images. IEEE Trans Med Imaging 18. Lee A, Schnabel J, Rajagopal V, Nielsen P, Nash M (2010) Breast 18(8):712–21 image registration by combining finite elements and free-form 25. Samani A, Bishop J, Yaffe MJ, Plewes DB (2001) Biomechanical deformations. In: Digital mammography, pp 736–743 3-D finite element modeling of the human breast using MRI data. 19. Liu GR, Quek SS (2013) The finite element method: a practical IEEE Trans Med Imaging 20(4):271–279 course. Butterworth-Heinemann, Boston 26. Schnabel JA, Tanner C, Castellano-Smith AD, Degenhard A, Leach 20. Maciel A, Boulic R, Thalmann D (2003) Deformable tissue MO, Hose DR, Hill DL, Hawkes DJ (2003) Validation of non- parameterized by properties of real biological tissue. In: Surgery rigid image registration using finite-element methods: application simulation and soft tissue modeling, pp 74–87. Springer to breast MR images. IEEE Trans Med Imaging 22(2):238–247 21. Pathmanathan P, Gavaghan DJ, Whiteley JP, Chapman SJ, Brady 27. Whelan B, Liney GP, Dowling JA, Rai R, Holloway L, McGarvie JM (2008) Predicting tumor location by modeling the deformation L, Feain I, Barton M, Berry M, Wilkins R, Keall P (2017) An MRI- of the breast. IEEE Trans Biomed Eng 55(10):2471–2480 compatible patient rotation system design, construction, and first 22. Picinbono G, Delingette H, Ayache N (2003) Non-linear organ deformation results. Med Phys 44(2):581–588 anisotropic elasticity for real-time surgery simulation. Graph Mod- els 65(5):305–321 23. Roose L, De Maerteleire W, Mollemans W, Suetens P (2005) Validation of different soft tissue simulation methods for breast augmentation. In: International congress series, vol 1281, pp 485– 490. Elsevier http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png International Journal of Computer Assisted Radiology and Surgery Springer Journals

Analytical derivation of elasticity in breast phantoms for deformation tracking

Download PDF
10 pages

Loading next page...
 
/lp/springer_journal/analytical-derivation-of-elasticity-in-breast-phantoms-for-deformation-sHerosOsNd
Publisher
Springer Journals
Copyright
Copyright © 2018 by The Author(s)
Subject
Medicine & Public Health; Imaging / Radiology; Surgery; Health Informatics; Computer Imaging, Vision, Pattern Recognition and Graphics; Computer Science, general
ISSN
1861-6410
eISSN
1861-6429
DOI
10.1007/s11548-018-1803-x
pmid
29869320
Publisher site
See Article on Publisher Site

Abstract

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. Keywords Biopsy · Magnetic resonance imaging · Elastic calibration · Breast Introduction different acquisition modalities and includes mammography (X-ray), ultrasound (US) and MRI. Screening and staging of breast cancer for diagnosis and sub- After image acquisition, proper localization of the tumor sequent treatment is based on medical images acquired on is essential for biopsy procedures to take tissue samples or to remove the tumor during surgery. To take full benefit from This project has received funding from the European Unions Horizon the previously acquired medical images, the location of the 2020 research and innovation programme under Grant Agreement No. tumor should be aligned from the preoperative imaging into the operating room. The position of the patient can vary from prone during MRI scanning to supine position required B Vincent Groenhuis v.groenhuis@utwente.nl for breast surgery for example. During ultrasound scanning and ultrasound-guided biopsy, the patient is returned on her Francesco Visentin francesco.visentin@univr.it back and additional compression is induced by the ultrasound probe. The computation of the elastic properties will serve University of Twente, Drienerlolaan 5, 7522 NB Enschede, as input for real-time adjustments of realistic deformations The Netherlands between preoperative and intra-operative images. For effec- University of Verona, Strada le Grazie 15, 37134 Verona, Italy 123 1642 International Journal of Computer Assisted Radiology and Surgery (2018) 13:1641–1650 tive deformation models, the elasticity of the model needs based image registration techniques, as in [9]or[18]. The to be known with good accuracy, i.e., the difference between sliding motion of the breast on the chest wall was observed computed and actual elasticity must be small. In this study, we [6], but usually a fixed muscle surface is applied during the aim for a maximum difference in the order of 10%, or at most FEM simulations [14,18,21]. two times the elasticity variation among FEM-simulated elas- This study introduces a method to analytically derive the ticity values. Image registration techniques based on image elastic modulus of the breast from a pair of MRI scans, taking intensities could be used for small deformations [24], but do local differences in tissue density and elasticity into account. not work in cases with large deformations such as the align- The two MRI scans differ by the direction of the gravitational ment from prone to supine configurations [4]. field, which are opposite to each other. Contrary to FEM- Deformation of the breast occurs due to body move- based numerical simulations, it is not needed to convert the ments. Various physics-based numerical procedures have MRI scan into a volumetric mesh, so mechanical properties been presented for biomechanical modeling and soft tissue on voxel scale are preserved. Also, only one iteration over all deformation. The most common computational schemes are voxels is necessary, which makes the method relatively fast based on linear or nonlinear biomechanical models including The proposed analytical method requires the breast to be mass-spring methods (MSM) [2,7,20,23], the mass-tensor vertically supported by a rigid planar base. As the rib cage method [10,22], the boundary element method [13,17] and is approximately cylindrical, a human breast would need to conventional finite element modeling (FEM) [3,25,26]. be supported by a patient-mounted flat plate with a hole for In an MSM system, an object is modeled by a collection the breast. In an MRI scanner, the breast coil could serve this of point masses linked together with massless springs. purpose. Recent studies show the use of FEM to align data with To avoid introduction of significant non-gravity-induced large deformations of the breast [15,16]. In FEM, a body is deformations when converting from prone to supine posi- subdivided into a set of finite elements (e.g., tetrahedral or tion, it is desirable to use a patient rotation system (PRS) hexahedra in 3D, triangles or other polygons in 2D). Dis- that allows leaving the patient on the bed with breast coil placements and positions of each element are approximated attached, while being flipped over by 180°. Such a system from discrete nodal values using interpolation functions: has been developed previously by Whelan et al. [27], which theoretically could be used to take MRI scans of a planar- supported breast in both prone and supine position. It may φ(x ) = h (x )φ (1) i i also be possible to tilt certain MRI scanners such as the 0.25 T G-scan Brio (Esaote SpA, Genoa, Italy), although this is where h is the interpolation function for the element con- generally limited to rotation over 90°only. taining x and φ is the scalar weight associated with h . i i Different choices for the element type and the interpolation functions exist, which depend on the accuracy requirements, Materials and methods geometry of the objects and computational complexity [19]. In general, FEM is used to solve a dynamic problem, which Four breast phantoms were constructed (Fig. 1, right), con- is expressed as partial differential equations (PDEs). These sisting of a rigid base with three fiducials, stiff superficial PDEs are then approximated with FEM. The FEM pro- tissue, soft deep tissue and 3–4 lesions. cedure has the advantage that it can handle complicated The superficial and deep tissues and lesions were made of geometries (and boundaries) of high quality. A dataset of polyvinyl chloride (PVC) with plasticizer mixed in different radiological 3D images of the breast anatomy (computed ratios to obtain different stiffnesses. Contrary to gelatin- tomography (CT) or MRI) is required to generate a patient- based phantoms, PVC is a durable material that can stay specific FEM. An advantage of MRI is that it shows high intact for extended periods. The superficial tissue consists of sensitivity for detecting breast tumors [8]. ThemainFEM relatively stiff PVC which was shaped using a pair of molds steps include: tissue classification/segmentation, tissue sur- (Fig. 1, left) and afterward filled with soft PVC to mimic deep face reconstruction, FEM volumetric mesh generation and tissue. The lesions were cut in different sizes and shapes from tissue type assignment for the FEM mesh. a block of relatively stiff PVC, placed inside the deep tissue at A patient-specific biomechanical model [11] was pre- random locations. A rigid frame was put on top and covered sented before to provide an initial deformation of the breast with a layer of stiff PVC. The four phantoms which were before registration between prone and supine MRI images. manufactured this way differ only in the stiffness of deep A zero-gravity reference state for both prone and supine tissue and the placement of lesions. configurations was estimated. The patient-specific unloaded Figure 2 shows the outline of a breast phantom in a neu- configuration was obtained [12]. The biomechanical meth- tral reference state. Depending on the orientation (prone or ods serve in most cases for the initialization of intensity- supine), it is deformed by the gravitational field and tip is 123 International Journal of Computer Assisted Radiology and Surgery (2018) 13:1641–1650 1643 Fig. 1 Left: pair of molds (yellow, green) for manufacturing superficial tissue (red). Right: one PVC breast phantom mounted in prone position 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. Each of the four phantoms was scanned in a 0.25 T MRI scanner (G-Scan Brio) using the 3D balanced steady- state free precession (bSSFP) sequence, with parameters TR=10ms, TE=5ms,FA = 60 , acquisition resolution 1.5 × 1.8 × 2.0 mm and isotropic reconstruction resolution 0.94 mm. Fig. 2 Breast in coil, with gravity-induced deformations in prone and supine positions (dashed lines) 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. displaced toward the anterior or posterior direction. The mag- The ideal, distorted and corrected grid patterns are shown in nitude of these deformations is related to the elasticity, and Fig. 3. The measured residual error is 0.2 mm, so sub-pixel the approach of the research is to reconstruct the elasticity resolution is feasible. from these deformations using different methods. Fig. 3 Left: MRI calibration grid. Right: actual (yellow), observed (blue) and distortion-corrected (red) grid locations of the calibration cube 123 1644 International Journal of Computer Assisted Radiology and Surgery (2018) 13:1641–1650 Fig. 4 Left: Example sagittal MRI slice. Right: Phantom I in prone and supine configuration, superimposed The distortion-corrected MRI scans (Fig. 4, left) were seg- mented 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. Elasticity estimation Preamble The deformation of an object in a gravitational field is the Fig. 5 Schematic view of force and pressure at a given height result of elongations of tissue, which depends on the local ratio of tensile stress σ and Young’s modulus E: section) that the deformation displacement can be solved analytically. The stress at a given location is primarily induced by the We introduce the assumption that the tensile stress σ weight of the masses below that location, and also influenced solely depends on the vertical position in the object, i.e., by interactions with surrounding tissue. In the general case, it is constant within any planar cross section parallel to the resulting stress distribution in the object is a complex the base. It can be shown that this assumption is valid for pattern and cannot be solved analytically, requiring simu- blocks, cylinders and prism-shaped objects which have a lations to quantify the deformations. However, in our case, constant cross section. For the breast phantom shapes, the we can use the knowledge that the object’s attachment to assumption can be justified by the fact that the masses the rigid frame is planar and perpendicular to the gravity of the whole breast are substantially positioned below the direction, when in prone and supine positions. For objects rigid base. To validate this assumption, the stress distribu- with a constant cross section such as a block or a cylinder, tion and elongation for a range of geometric shapes are also it can be shown (see “Analytical derivation of elasticity” investigated. 123 International Journal of Computer Assisted Radiology and Surgery (2018) 13:1641–1650 1645 Analytical derivation of elasticity The displacement equation can now be written as follows: g m(h) Figure 5 schematically shows the forces and pressures acting D = dh (8) on a shape with inhomogeneous density and elasticity, hang- E ˆ 0 A(h)E (h) ing from a planar, rigid attachment on the top. At a given height h, the cross-sectional area is A(h), the mass of the It can split into an object-specific intrinsic part which body below it is denoted as m(h) and the gravitational force remains constant across all simulations, and an extrinsic (variable) part depending on g and E only. The intrinsic part acting on it F (h). We now derive expressions for the vertical stress σ(h) and elongation (h) for every height, leading to β is defined as: a formula for the displacement D of the lower extremity of m(h) the body. β = dh (9) The total mass of the body up to height h is given as: ˆ 0 A(h)E (h) h Substituting into D gives: m(h) = ρ(x , y, z)dxdydz (2) D = β (10) The gravitational force acting on the slice at height h is For the scanned breast phantoms, we, therefore, assume calculated as: that the displacement (for small displacements) is linear in g/E, with proportionality factor β.The β value can be esti- F (h) = m(h)g (3) mated from DICOM data, in combination with knowledge of the materials. For PVC phantoms, its density was measured The tensile stress in the slice is generally not constant, to be ρ = 1.075 g/cm . and its exact distribution depends on many levels of tissue Analyzing the prone and supine scans of a phantom, we interactions. We are interested in the mean tensile stress σ(h), have β and β for prone and supine, respectively. In gen- p s which is found by dividing the gravitational force by the eral, β = β , because the shapes are significantly different: p s slice’s cross-sectional area: The total volume and cross-sectional area at the base are approximately equal, but due to difference in height the cross- F (h) m(h)g sectional shape is more squeezed in prone position than in σ(h) = = (4) A(h) A(h) the supine one. The phantom height H is ill-defined due to possible irreg- The tissue elasticity is also inhomogeneous in general, ularities at the tip, but the difference H = H − H can p s with local Young’s modulus E (r), again averaged to E (h) be accurately determined by comparing point clouds around for height h. The local relative elongation is  = L/L = the tip using, e.g., the iterative closest point algorithm [5], σ(r)/E (r), and the mean elongation at height h is given as: and optimizing H such that the total point distance is min- imal, or alternatively by comparing the centroids of the point σ(h) m(h)g clouds. (h) = = (5) E (h) A(h)E (h) The parameter we want to compute is the Young’s mod- ulus E. When no forces act on the phantom, it would have The total displacement of the body’s lower extremity is some shape halfway the prone and supine shapes. The tip dis- found by integrating all infinitesimal elongations: placement to either prone or supine shape in a gravitational field g,is H /2. We can now derive the Young’s modulus H H E as follows: m(h) D = (h)dh = g dh (6) A(h)E (h) 0 0 β + β p s β = (11) The purpose of this study is to find the average Young’s β g 2(β + β )g n p s E = = (12) modulus E from a pair of gravity-induced body displace- H /2 H ments. To preserve differences in (mean) elasticity among slices, we factorize every slice’s elasticity into a constant Numerical simulation of deformations factor E and a layer-specific adjustment factor E (h): The purpose of FEM simulations is to determine the elas- E (h) = E E (h) (7) ticity E of the different phantoms, based on the segmented 123 1646 International Journal of Computer Assisted Radiology and Surgery (2018) 13:1641–1650 models. The general strategy is to apply a gravitational field value. The mean value (square-harmonic mean-root) of E sp to the FEM model of a phantom in a specific direction. This and E is then taken as the elasticity of the final phantom. ps deformed model is then compared to a reference phantom which was scanned in a different orientation, providing infor- mation about the elasticity parameter. Results In the following subsections, we present two strategies to find the Young’s modulus by simulation, of which one Validation of analytical stress calculation on strategy is performed by two different simulation software geometric shapes packages. Nine homogeneous geometric shapes were generated and Estimating the ˇ values by simulation in SOFA analyzed: two cylinders with different aspect ratios, a cone, a T-piece in normal and upside-down orientation, a half sphere, In “Analytical derivation of elasticity” section, we have a sphere, an hourglass and a snake-like shape. introduced a method to derive the values of β for the four Figure 6 shows the stress distribution along the vertical phantoms in different orientations directly from a DICOM midway plane for all nine shapes. The first row uses the ana- scan. In this section, we find β by simulation in SOFA at five lytical computation method. The assumption that the stress different mesh resolutions [1]. For each mesh resolution, we distribution is constant in a cross-sectional area parallel to the have run a simulation with the phantom’s Young’s modulus base, is reflected in having constant colors in horizontal direc- set to E = 6000 Pa and gravity g = 2.0m/s . After 100 iter- tion. The second row shows the tensile stress from numerical ations, the simulation has stabilized and the vertices of the simulations using the SOFA software package under the same mesh in this configuration were extracted and analyzed. The conditions. displacement from the initial position follows by comparing Table 1 lists the calculated and simulated β values for the the point clouds around the tip. The value of β then follows same geometric shapes. from Eq. (10). This procedure is repeated for each resolution The following observations can be made: of the mesh and for both prone and supine orientations, then the mean β and β values were computed. From the β , β s p s p – For cylinder, cubic and prism-like shapes that have a and H , and assuming linearity of the displacement to g/E constant cross-sectional area (a and b), the numeri- ratio, the Young’s modulus E can be derived using Eqs. (11) cally derived stress distribution matches the analytically and (12). derived one quite well. The β values derived by both methods are well comparable (deviation under 9%). Supine–prone and prone–supine simulation and matching – For shapes that do not have a constant cross-sectional in SOFA and Febio area, but are substantially vertically supportive (c-h), the analytically calculated and SOFA-simulated β values are Taking a phantom scanned in supine configuration, the base still comparable (deviation up to 18%) although the stress of the phantom is immobilized and a force field sized two distribution is different. times the gravity (19.62 m/s ) in anterior direction is applied – For shapes in which the lower extremity is not vertically to the phantom. After stabilization in simulation, the final supported by the base, i.e., no vertical line of maximum state is extracted and compared to the phantom in prone posi- height can be drawn that entirely lies within the model tion, which serves as the reference phantom. (i), both the analytically calculated β value and the stress The error value, , is defined as the distance between the distribution are inconsistent with simulations. 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. Analytical derivation of elasticity of phantoms The minimization is performed using the Newton’s method computed over E and the distance error, corrected Each of the four phantoms was scanned in prone and supine by an adaptive step approach (when the FEM analysis soft- position, and from the resulting DICOM scans, the β and ware diverges). When procedure ends, i.e., when the method β values are computed using Eq. 9 and assuming a homoge- achieves a pre-defined error or when it reaches a maximum neous density and elasticity distribution. From these values number of iterations, the estimated E parameter is returned plus the observed vertical displacements, the E parameters with its associated error. are computed using Eq. 12 and the results are listed in Table 2. The procedure is then repeated for the opposite direction It can be observed that phantom IV has the highest β and E (prone to supine). In general, this also leads to a different E values, making it the stiffest phantom, while phantom II is 123 International Journal of Computer Assisted Radiology and Surgery (2018) 13:1641–1650 1647 Fig. 6 Analytically derived tensile stress (top row) compared with simulated stress (bottom row) for a selection of geometric shapes Table 1 Calculated and simulated β values for the nine geometric Table 3 Properties of four phantoms, derived by numerical simulation shapes in SOFA in five different resolution scales and then averaged Geometric shape Calculated β Simulated β Phantom β β HE s p a 2375 2169 I 1007 ± 58 1134 ± 38 3.28 6403 ± 207 b 2373 2229 II 947 ± 41 1125 ± 36 4.73 4297 ± 113 c 772 724 III 1131 ± 48 1259 ± 61 3.58 6549 ± 213 d 1638 1581 IV 1170 ± 43 1383 ± 34 2.93 8548 ± 184 e 4500 4979 f 213 205 g 1932 2276 Figure 7 shows the analytically derived stress distribution h 3802 3797 in the transversal plane of phantom I in supine configuration i 4942 26,499 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 cer- Table 2 Analytically derived properties of four phantoms, under the tain geometric shapes in Fig. 6a–h. Only the analytic method assumption of constant tensile stress in each cross section shows a sharp transition at the boundary layer, as the ana- Phantom β β HE s p lytical method uses slices with thickness of one voxel while the FEM-based method subdivides the volume in a different I 1215 1298 3.28 7514 way. II 1129 1269 4.73 4972 III 1356 1444 3.58 7673 IV 1420 1471 2.93 9677 Numerical simulation by supine–prone and prone–supine matching in SOFA Table 4 lists the elasticities obtained by numerical simulation the softest one. In general, the β values are higher in prone from supine to prone position and vice versa, in SOFA. As position, which is as expected. opposed to the β computation method, the prone–supine sim- ulation method also takes nonlinearities into account which Simulation of ˇ in SOFA theoretically results in a more accurate estimate of the E value. For numerical FEM simulations, each DICOM scan was For each resolution, up to ten simulation runs are needed segmented and meshed at five different levels of detail and to find the final E value in which the error vanishes. This subsequently simulated in the SOFA simulation package. The makes the method relatively slow, requiring about twenty resulting β values of the four phantoms (in both orientations) minutes of computation time on a quad-core 2.5 GHz com- plus the averaged E value are listed in Table 3. Calculation puter per phantom. By parallelizing computations of the four of each β value requires ten simulation runs in SOFA, lasting phantoms, the total computation time for all E values was a few minutes in total. measured to be approximately half an hour. 123 1648 International Journal of Computer Assisted Radiology and Surgery (2018) 13:1641–1650 Fig. 7 Tensile stress for phantom I in the transversal plane, in supine position. Left: derived using analytical method. Center and right: numerically simulated using SOFA in low resolution (center) and high resolution (right). The dashed line indicates the boundary plane between the rigid and deformable parts Table 4 Elasticity values found by numerical simulations from supine- to-prone (E ) and prone-to-supine (E ) in four different resolution sp ps scales and then averaged, using SOFA Phantom E E Mean E sp ps I 5047 ± 374 6459 ± 373 5688 ± 272 II 3395 ± 189 4513 ± 272 3895 ± 159 III 5381 ± 376 6828 ± 438 6040 ± 288 IV 6245 ± 433 7445 ± 322 6805 ± 283 Table 5 Elasticity values found by simulating from supine-to-prone (E ) and prone-to-supine (E ) in four different resolution scales and sp ps then averaged, using FEBio as software package Phantom E E Mean E Fig. 8 Young’s modulus for four phantoms, derived by four different sp ps methods I 5046 ± 272 5252 ± 307 5145 ± 254 II 4290 ± 351 4298 ± 273 4291 ± 276 III 5291.52 ± 383 5639 ± 456 5459 ± 400 Table 6 Mean elasticity values Phantom Mean E for each phantom, taken as the IV 7916 ± 1165 7564 ± 957 7731 ± 1016 average of the separate values I 6188 ± 886 derived by the four different II 4364 ± 387 methods III 6430 ± 815 IV 8190 ± 1057 Numerical simulation by supine–prone and prone–supine matching in FEBio Table 5 lists the elasticity values using the FEBio software package. The resulting elasticity values are comparable to The second method involves finding E directly by simulation those obtained by SOFA. A relatively high variance is present from supine to prone position such that the tip position error in phantom IV, which may be caused by side effects in the is eliminated. The first method seems to give consistently software package. higher estimates for E, especially for phantom IV. Possible causes might be the nonlinearity of the displacement-to- Comparison of different elasticity measurement g/E ratio, i.e., β cannot be considered constant for the methods required range of displacements. Furthermore, the defor- mations of the tip resulting from proper FEM simulations Figure 8 graphically shows the elasticity of the four phan- influence the displacement calculations. As the second algo- toms, derived using the different methods, while Table 6 lists rithm uses the iterative point cloud algorithm to minimize tip the overall phantom elasticities, averaged from the four dif- displacements and also takes nonlinear effects into account, ferent methods. that one can be considered more accurate than the first Two methods using SOFA were presented: The first one one. numerically simulates the β and β values from the supine The numerical results from FEBio simulations are in s p and prone meshes separately and measures the tip displace- accordance with SOFA matching simulations, which is an ment D, from which the phantom’s elasticity E is derived. indication that the simulations are consistent. 123 International Journal of Computer Assisted Radiology and Surgery (2018) 13:1641–1650 1649 Discussion References 1. Allard J, Cotin S, Faure F, Bensoussan PJ, Poyer F, Duriez C, We have presented a new method to analytically evaluate the Delingette H, Grisoni L (2007) Sofa-an open source framework for elasticity of breast phantoms, from a pair of MRI scans in medical simulation. In: MMVR 15-medicine meets virtual reality, prone and supine position. The values found from analyz- vol 125, pp 13–18. IOP Press ing the gravity-induced deformations are comparable to the 2. Altomonte M, Zerbato D, Botturi D, Fiorini P (2008) Simula- tion of deformable environment with haptic feedback on GPU. elasticities derived from FEM simulations using FEBio and In: IEEE/RSJ international conference on intelligent robots and SOFA, with deviations of up to 18%. A study on nine geomet- systems, 2008. IROS 2008, pp 3959–3964. IEEE ric shapes has shown that the method is not only applicable to 3. Azar FS, Metaxas DN, Schnall MD (2001) A deformable finite breast shapes, but also to other bodies and geometric objects element model of the breast for predicting mechanical deformations under external perturbations. Acad Radiol 8(10):965–975 as long as it is substantially supported by a planar rigid base. 4. Behrenbruch C, Marias K, Armitage P, Moore N, Clarke J, Brady The advantages of the analytical method are that the J (2001) Prone-supine breast MRI registration for surgical visual- elasticity calculation is very fast (< 1 s) and takes each isation. In: Medical image understanding and analysis individually scanned voxel into account, without need for 5. Besl PJ, McKay ND (1992) A method for registration of 3-D shapes. IEEE Trans Pattern Anal Mach Intell 14(2):239–256 mesh generation. As the voxel intensity in a scan gives cer- 6. Carter T, Tanner C, Beechey-Newman N, Barratt D, Hawkes D tain information about tissue type, density and/or elasticity (2008) MR navigated breast surgery: method and initial clinical (depending on scanning protocol), tissue inhomogeneities experience. In: Medical image computing and computer-assisted can be directly incorporated in the analytical computations. intervention-MICCAI 2008, pp 356–363 7. Chang YH, Chen YT, Chang CW, Lin CL (2010) Development The main limitations are that the method is only suitable for scheme of haptic-based system for interactive deformable simula- deformations in the linear range, and that the shapes must be tion. Comput Aided Des 42(5):414–424 substantially supported by a planar base perpendicular to the 8. Chevrier MC, David J, Khoury ME, Lalonde L, Labelle M, Trop gravitational field. I (2016) Breast biopsies under magnetic resonance imaging guid- ance: challenges of an essential but imperfect technique. Curr Probl The fact that a human breast is relatively flexible and the Diagn Radiol 45(3):193–204. https://doi.org/10.1067/j.cpradiol. chest wall is not planar but cylindrically shaped, makes clin- 2015.07.002 ical application difficult. An artificial planar support base 9. Conley RH, Meszoely IM, Weis JA, Pheiffer TS, Arlinghaus LR, could be constructed by using a patient-mounted breast coil, Yankeelov TE, Miga MI (2015) Realization of a biomechanical model-assisted image guidance system for breast cancer surgery ideally in combination with a patient rotation system. The using supine MRI. Int J Comput Assist Radiol Surg 10(12):1985– presented methods may also have applications in different domains, wherever deformation of bodies is involved in sit- 10. Cotin S, Delingette H, Ayache N (2000) A hybrid elastic model uations that meet the aforementioned boundary conditions for real-time cutting, deformations, and force feedback for surgery training and simulation. Vis Comput 16(8):437–452 The conclusion is that under specific conditions, the elas- 11. Eiben B, Han L, Hipwell J, Mertzanidou T, Kabus S, Bülow T, ticity of a deformable object such as a human breast can Lorenz C, Newstead G, Abe H, Keshtgar M, Ourselin S, Hawkes DJ be quickly computed from a pair of volumetric scans with (2013) Biomechanically guided prone-to-supine image registration sufficient accuracy, without need for FEM simulations. This of breast MRI using an estimated reference state. In: 2013 IEEE 10th international symposium on biomedical imaging (ISBI), pp promising result opens the door to new applications which 214–217. IEEE can benefit from this complementary and near-real-time elas- 12. Eiben B, Vavourakis V, Hipwell JH, Kabus S, Lorenz C, Buelow ticity computation method. T, Hawkes DJ (2014) Breast deformation modeling: comparison of methods to obtain a patient specific unloaded configuration. In: Proceedings of SPIE, vol 9036, pp 903615–903618 Compliance with ethical standards 13. Greminger MA, Nelson BJ (2003) Deformable object tracking using the boundary element method. In: 2003 IEEE computer soci- ety conference on computer vision and pattern recognition, 2003. Proceedings, vol 1, pp I–I. IEEE Conflict of interest The authors declare that they have no conflict of 14. Han L, Hipwell J, Mertzanidou T, Carter T, Modat M, Ourselin S, interest. Hawkes D (2011) A hybrid fem-based method for aligning prone and supine images for image guided breast surgery. In: 2011 IEEE Human and animal rights This article does not contain any studies with international symposium on biomedical imaging: from nano to human participants or animals performed by any of the authors. This macro, pp 1239–1242. IEEE articles does not contain patient data. 15. Han L, Hipwell JH, Eiben B, Barratt D, Modat M, Ourselin S, Hawkes DJ (2014) A nonlinear biomechanical model based regis- Open Access This article is distributed under the terms of the Creative tration method for aligning prone and supine MR breast images. Commons Attribution 4.0 International License (http://creativecomm IEEE Trans Med Imaging 33(3):682–694 ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, 16. Han L, Hipwell JH, Tanner C, Taylor Z, Mertzanidou T, Cardoso J, and reproduction in any medium, provided you give appropriate credit Ourselin S, Hawkes DJ (2011) Development of patient-specific to the original author(s) and the source, provide a link to the Creative biomechanical models for predicting large breast deformation. Commons license, and indicate if changes were made. Phys Med Biol 57(2):455 123 1650 International Journal of Computer Assisted Radiology and Surgery (2018) 13:1641–1650 17. James DL, Pai DK (2005) A unified treatment of elastostatic con- 24. Rueckert D, Sonoda LI, Hayes C, Hill DL, Leach MO, Hawkes tact simulation for real time haptics. In: ACM SIGGRAPH 2005 DJ (1999) Nonrigid registration using free-form deformations: courses, p 141. ACM application to breast MR images. IEEE Trans Med Imaging 18. Lee A, Schnabel J, Rajagopal V, Nielsen P, Nash M (2010) Breast 18(8):712–21 image registration by combining finite elements and free-form 25. Samani A, Bishop J, Yaffe MJ, Plewes DB (2001) Biomechanical deformations. In: Digital mammography, pp 736–743 3-D finite element modeling of the human breast using MRI data. 19. Liu GR, Quek SS (2013) The finite element method: a practical IEEE Trans Med Imaging 20(4):271–279 course. Butterworth-Heinemann, Boston 26. Schnabel JA, Tanner C, Castellano-Smith AD, Degenhard A, Leach 20. Maciel A, Boulic R, Thalmann D (2003) Deformable tissue MO, Hose DR, Hill DL, Hawkes DJ (2003) Validation of non- parameterized by properties of real biological tissue. In: Surgery rigid image registration using finite-element methods: application simulation and soft tissue modeling, pp 74–87. Springer to breast MR images. IEEE Trans Med Imaging 22(2):238–247 21. Pathmanathan P, Gavaghan DJ, Whiteley JP, Chapman SJ, Brady 27. Whelan B, Liney GP, Dowling JA, Rai R, Holloway L, McGarvie JM (2008) Predicting tumor location by modeling the deformation L, Feain I, Barton M, Berry M, Wilkins R, Keall P (2017) An MRI- of the breast. IEEE Trans Biomed Eng 55(10):2471–2480 compatible patient rotation system design, construction, and first 22. Picinbono G, Delingette H, Ayache N (2003) Non-linear organ deformation results. Med Phys 44(2):581–588 anisotropic elasticity for real-time surgery simulation. Graph Mod- els 65(5):305–321 23. Roose L, De Maerteleire W, Mollemans W, Suetens P (2005) Validation of different soft tissue simulation methods for breast augmentation. In: International congress series, vol 1281, pp 485– 490. Elsevier

Journal

International Journal of Computer Assisted Radiology and Surgery – Springer Journals

Published: Jun 4, 2018

References

  • A deformable finite element model of the breast for predicting mechanical deformations under external perturbations
    Azar, FS; Metaxas, DN; Schnall, MD
  • A method for registration of 3-D shapes
    Besl, PJ; McKay, ND
  • Medical Image Computing and Computer-Assisted Intervention – MICCAI 2008
    Carter, Timothy; Tanner, Christine; Beechey-Newman, Nicolas; Barratt, Dean; Hawkes, David
  • Development scheme of haptic-based system for interactive deformable simulation
    Chang, YH; Chen, YT; Chang, CW; Lin, CL
  • Breast biopsies under magnetic resonance imaging guidance: challenges of an essential but imperfect technique
    Chevrier, MC; David, J; Khoury, ME; Lalonde, L; Labelle, M; Trop, I
  • Realization of a biomechanical model-assisted image guidance system for breast cancer surgery using supine MRI
    Conley, RH; Meszoely, IM; Weis, JA; Pheiffer, TS; Arlinghaus, LR; Yankeelov, TE; Miga, MI
  • A hybrid elastic model for real-time cutting, deformations, and force feedback for surgery training and simulation
    Cotin, S; Delingette, H; Ayache, N
  • A nonlinear biomechanical model based registration method for aligning prone and supine MR breast images
    Han, L; Hipwell, JH; Eiben, B; Barratt, D; Modat, M; Ourselin, S; Hawkes, DJ
  • Development of patient-specific biomechanical models for predicting large breast deformation
    Han, L; Hipwell, JH; Tanner, C; Taylor, Z; Mertzanidou, T; Cardoso, J; Ourselin, S; Hawkes, DJ
  • Digital Mammography
    Lee, Angela W. C.; Schnabel, Julia A.; Rajagopal, Vijay; Nielsen, Poul M. F.; Nash, Martyn P.
  • The finite element method: a practical course
    Liu, GR; Quek, SS
  • Predicting tumor location by modeling the deformation of the breast
    Pathmanathan, P; Gavaghan, DJ; Whiteley, JP; Chapman, SJ; Brady, JM
  • Non-linear anisotropic elasticity for real-time surgery simulation
    Picinbono, G; Delingette, H; Ayache, N
  • Nonrigid registration using free-form deformations: application to breast MR images
    Rueckert, D; Sonoda, LI; Hayes, C; Hill, DL; Leach, MO; Hawkes, DJ
  • Biomechanical 3-D finite element modeling of the human breast using MRI data
    Samani, A; Bishop, J; Yaffe, MJ; Plewes, DB
  • Validation of nonrigid image registration using finite-element methods: application to breast MR images
    Schnabel, JA; Tanner, C; Castellano-Smith, AD; Degenhard, A; Leach, MO; Hose, DR; Hill, DL; Hawkes, DJ
  • An MRI-compatible patient rotation system design, construction, and first organ deformation results
    Whelan, B; Liney, GP; Dowling, JA; Rai, R; Holloway, L; McGarvie, L; Feain, I; Barton, M; Berry, M; Wilkins, R; Keall, P