Modeling of Personalized Anatomy using Plastic Strains
MModeling of Personalized Anatomy using Plastic Strains
BOHAN WANG,
University of Southern California, USA
GEORGE MATCUK,
University of Southern California, USA
JERNEJ BARBIČ,
University of Southern California, USA
Template (b) MRI slice
Our output
MRI sliceICP markers (d) maximum principal stretch (c) our output shape(a) template shape (e) principal stretch ratio (anisotropicity)
Fig. 1.
We optimized the shape of this hand muscle (Adductor Pollicis) to match the MRI scan. (a) template shape [C. Erolin 2019]; (b) representativeMRI slice [Wang et al. 2019]; we rigidly aligned the template mesh onto the markers in the MRI scan, producing the white contour which is obviously incorrect(deeply penetrates the bone; and extends out of the volume of the MRI-scanned hand); (c) our output shape, optimized to the MRI scan; the white contour nowmatches the scan; (d,e) large anisotropic spatially varying strains accommodated by our method (d: maximum principal stretch; e: ratio between maximumand minimum principal stretch).
We give a method for modeling solid objects undergoing large spatiallyvarying and/or anisotropic strains, and use it to reconstruct human anatomyfrom medical images. Our novel shape deformation method uses plasticstrains and the Finite Element Method to successfully model shapes under-going large and/or anisotropic strains, specified by sparse point constraintson the boundary of the object. We extensively compare our method to stan-dard second-order shape deformation methods, variational methods andsurface-based methods and demonstrate that our method avoids the spiki-ness, wiggliness and other artefacts of previous methods. We demonstratehow to perform such shape deformation both for attached and un-attached(“free flying”) objects, using a novel method to solve linear systems withsingular matrices with a known nullspace. While our method is applicableto general large-strain shape deformation modeling, we use it to create per-sonalized 3D triangle and volumetric meshes of human organs, based on
Authors’ addresses: Bohan Wang, University of Southern California, Los Angeles, CA,USA, [email protected]; George Matcuk, University of Southern California, LosAngeles, CA, USA, [email protected]; Jernej Barbič, University of Southern California,Los Angeles, CA, USA, [email protected] to make digital or hard copies of part or all of this work for personal orclassroom use is granted without fee provided that copies are not made or distributedfor profit or commercial advantage and that copies bear this notice and the full citationon the first page. Copyrights for third-party components of this work must be honored.For all other uses, contact the owner/author(s).© 2020 Copyright held by the owner/author(s).XXXX-XXXX/2020/8-ARThttps://doi.org/10.1145/nnnnnnn.nnnnnnn
MRI or CT scans. Given a medically accurate anatomy template of a genericindividual, we optimize the geometry of the organ to match the MRI or CTscan of a specific individual. Our examples include human hand muscles, aliver, a hip bone, and a gluteus medius muscle (“hip abductor”).CCS Concepts: •
Computing methodologies → Shape modeling ; Phys-ical simulation .Additional Key Words and Phrases: shape deformation, large strain, FEM,plastic strain, optimization, MRI, CT, anatomy
ACM Reference Format:
Bohan Wang, George Matcuk, and Jernej Barbič. 2020. Modeling of Per-sonalized Anatomy using Plastic Strains. 1, 1 (August 2020), 18 pages.https://doi.org/10.1145/nnnnnnn.nnnnnnn
Modeling and simulating human anatomy is very important in manyapplications in computer graphics, animation, medicine, film andreal-time systems such as games and virtual reality. In this paper,we demonstrate how to model anatomically realistic personalized three-dimensional shapes of human organs, based on medical im-ages of a real person, such as Magnetic Resonance Imaging (MRI)or Computed Tomography (CT). Such modeling is crucially impor-tant for personalized medicine. For example, after scanning thepatient with an MRI or CT scanner, doctors can use the resulting , Vol. 1, No. 1, Article . Publication date: August 2020. a r X i v : . [ c s . G R ] A ug • Bohan Wang, George Matcuk, and Jernej Barbič
3D meshes to perform pre-operative surgery planning. Such modelsare also a starting point for anatomically based human simulationfor applications in computer graphics, animation and virtual reality.Constructing volumetric meshes that match an organ in a medicalimage can also help with building volumetric correspondences be-tween multiple MRI or CT scans of the same person [Rhee et al.2011], e.g., for medical education purposes.Although the types, number and function of organs in the hu-man body are generally the same for any human, the shape ofeach individual organ varies greatly from person to person, dueto the natural variation across the human population. The shapevariation is substantial: any two individuals’ organs Ω ⊂ R and Ω ⊂ R generally vary by a non-trivial shape deformation func-tion Φ : Ω → Ω that often consists of large and spatially varyinganisotropic strains (see Figure 1). By “large anisotropic strain”, wemean that the singular values of the 3x3 gradient matrix of Φ areboth different to each other and substantially different from 1.0, i.e.,the material locally stretches (or compresses) by large amounts; andthis amount is different in different directions and varies spatiallyacross the model.We tackle the problem of how to model such large shape vari-ations, using volumetric 3D medical imaging (such as MRI or CTscan), and a new shape deformation method capable of modelinglarge spatially varying anisotropic strains. We note that the bound-ary between the different organs in medical images is often blurry.For example, in an MRI of a human hand, the muscles often “blend”into each other and into fat without clear boundaries; a CT scan haseven less contrast. We therefore manually select as many reliablepoints ("markers") as possible on the boundary of the organ in themedical image; some with correspondence (“landmark constraints”)to the same anatomical landmark in the template organ, and somewithout (“ICP constraints”). Given a template volumetric mesh ofan organ of a generic individual, a medical image of the same organof a new individual, and a set of landmark and ICP (Iterative ClosestPoint) constraints, our paper asks how to deform the template meshto match the medical image.Our first attempt to solve this problem was to use standard shapedeformation methods commonly used in computer graphics, suchas as-rigid-as-possible energy (ARAP) [Sorkine and Alexa 2007],bounded biharmonic weights (BBW) [Jacobson et al. 2011], bihar-monic weights with linear precision (LBW) [Wang et al. 2015a], anda finite element method static solver (FEM) [Barbič et al. 2009]. Aswe show in Section 3.8, none of these methods was able to capturethe large strains observed in medical images. Namely, these stan-dard methods either cannot model point constraints when the shapeundergoes large spatially varying strains, or introduce excessivecurvature. For example, in the limit where the tetrahedral meshis refined to finer and finer tets, the FEM static solver produces aspike output (Figure 2); and similar limitations apply to the othermethods.Instead, we give a new shape deformation method that uses plas-tic strains and the Finite Element Method to successfully modelshapes undergoing large and/or anisotropic strains, controlled bythe sparse landmark and ICP point constraints on the boundary ofthe object. In order to do so, we formulate a nonlinear optimiza-tion problem for the unknown plastic deformation gradients of the Fig. 2.
Second-order methods produce spiky outputs as the tet meshis refined.
Here, we show the output of an FEM static solver under a singlehard point constraint (seen in (a)), i.e., we minimize the FEM elastic energyunder the shown hard point constraint. Bunny is fixed at the bottom. Pois-sons’s ratio is 0.49. As we increase the tet mesh resolution in (b)-(e), thespike becomes progressively narrower, which is undesirable. Changing theelastic stiffness (Young’s modulus) of the bunny changes nothing. Convert-ing the constraint into a soft spring constraint also does not help (f); now,the constraint is not even satisfied. template shape Ω , such that under these gradients, the shape Ω transforms into a shape that matches the medical image landmarkand ICP constraints. The ICP constraints are handled by properlyincorporating the ICP algorithm into our method. We note that insolid mechanics, plastic deformation gradients are a natural tool tomodel large volumetric shape variations. We are, however, unawareof any prior work that has used plastic deformation gradients andthe Finite Element Method to model large-strain shape deformation.In order to make our method work, we needed to overcome severalnumerical obstacles. The large-strain shape optimization problemis highly nonlinear and cannot be reliably solved with off-the-shelfoptimizers such as the interior point method [Artelys 2019]. Further-more, a naive solution requires solving large dense linear systemsof equations. We demonstrate how to adapt the Gauss-Newtonoptimization method to robustly and efficiently solve our shapedeformation problem, and how to numerically avoid dense linearsystems. In order to optimize our shapes, we needed to derive ana-lytical gradients of internal elastic forces and the tangent stiffnessmatrix with respect to the plastic strain, which will be useful forfurther work on using plasticity for optimization and design of 3Dobjects. In addition, we address objects that are attached to otherobjects, such as a hand muscle attached to one or more bones; aswell as un-attached objects. An example of an un-attached object isa liver, where the attachments to the surrounding tissue certainlyexist, but are not easy to define. It is practically easier to just modelthe liver as an un-attached object. In order to address un-attachedobjects, we give a novel numerical method to solve linear systemswith singular matrices with a known nullspace. Such linear sys-tems are commonly encountered in applications in geometric shape , Vol. 1, No. 1, Article . Publication date: August 2020. odeling of Personalized Anatomy using Plastic Strains • 3 modeling and nonlinear elastic simulation. Our examples includehuman hand muscles, a liver, a hip bone and a hip abductor muscle(“gluteus medius”), all of which undergo substantial and non-trivialshape change between the template and the medical image. In this section, we introduce closely related work and discuss therelationship to our work.
Geometric shape modeling.
Geometric shape modeling is an im-portant topic in computer graphics research; e.g., see the Botschand Sorkine [Botsch and Sorkine 2008] survey and the SIGGRAPHcourse notes by Alexa et al. [Alexa et al. 2006]. Popular methodsinclude variational methods [Botsch and Kobbelt 2004], Laplaciansurface editing [Sorkine et al. 2004], as-rigid-as-possible (ARAP)deformation [Igarashi et al. 2005; Sorkine and Alexa 2007], coupledprisms [Botsch et al. 2006] and partition-of-unity methods such asbounded biharmonic weights (BBW) [Jacobson et al. 2011] and bihar-monic weights with linear precision [Wang et al. 2015a]; we providea comparison in Section 3.8 and in several other Figures in the paper.Our method reconstructs the surface shape from a set of un-orientedpoint observations; this goal is similar to variational implicit surfacemethods [Huang et al. 2019; Turk and O’Brien 1999]; we give a com-parison in Section 4. Point clouds can also be used to optimize restshapes [Twigg and Kačić-Alesić 2011] and material properties of 3Dsolids [Wang et al. 2015b]. Such a method cannot be applied to ourproblem because it assumes a 4D dense point cloud input; whereaswe assume 3D sparse point inputs as commonly encountered inmedical imaging. Point constraint artifacts of second-order meth-ods can be addressed using spatial averaging [Bergou et al. 2007;Kavan et al. 2011]; however this requires specifying the averagingfunctions (often by hand) and, by the nature of averaging, causesthe constraints to be met only approximately. Our method can meetthe constraints very closely (under 0 . Plasticity.
Elastoplastic simulations are widely used in computeranimation. O’Brien et al. [2002] and Muller and Gross [2004] usedan additive plasticity formulation, whereas Irving et al. [Irving et al.2004] presented a multiplicative formulation and argued that it isbetter for handling large plasticity; we adopt multiplicative formu-lation in our work. The multiplicative model was used in manysubsequent publications to simulate plasticity, e.g., [Bargteil et al.2007; Chen et al. 2018; Stomakhin et al. 2013]. Because plasticitymodels shapes that undergo permanent and large deformation, itis in principle a natural choice also for geometric shape modeling.However, such an application is not straightforward: an incorrectchoice of the optimization energy will produce degenerate outputs,elastoplastic simulations in equilibrium lead to linear systems withsingular matrices, optimization requires second-order derivativesfor fast convergence, and easily produces large linear systems withdense matrices. We present a solution to these obstacles. To the bestof our knowledge, we are the first paper to present such a compre-hensive approach for using plasticity for geometric shape modelingwith large and anisotropic strains.
Anatomically based simulation.
Anatomically based simulation ofthe human body has been explored in multiple publications. Forexample, researchers simulated human facial muscles [Sifakis et al.2005], the entire upper human body [Lee et al. 2009], volumetricmuscles for motion control [Lee et al. 2018] and hand bones andsoft tissue [Wang et al. 2019]. Anatomically based simulation isalso popular in film industry [Tissue 2013]. Existing papers largelysimulate generic humans because it is not easy to create accurateanatomy personalized to each specific person. Our method canprovide such an input anatomy, based on a medical image of anyspecific new individual.
Medical image registration.
Deformable models are widely used inmedical image analysis [McInerney and Terzopoulos 2008]. Extract-ing quality anatomy geometry from medical images is difficult. Forexample, Sifakis and Fedkiw [2005] reported that it took them “sixmonths” (including implementing the tools) to extract the facialmuscles from the visible human dataset [U.S. National Library ofMedicine 1994], and even with the tools implemented it would stilltake “two weeks”. With our tools, we are able to extract all the 17muscles of the human hand in 1 day (including computer and user-interaction time). Bones generally have good contrast against thesurrounding tissue and can be segmented using active contour meth-ods [Székely et al. 1996] or Laplacian-based segmentation [Grady2006; Wang et al. 2019]. For bones, it is therefore generally possibleto obtain a “dense” set of boundary points in the medical image.Gilles et al. [Gilles et al. 2010] used this to deform template skele-ton models to match a subject-specific MRI scan and posture. Theyused the ARAP energy and deformed surface meshes. In contrast,we give a method that is suitable for soft tissues where the imagecontrast is often low (our hand muscles and liver examples) andthat accommodates volumetric meshes and large volumetric scalingvariations between the template and the subject. If one assumesthat the template mesh comes with a registered MRI scan (or ifone manually creates a template mesh that matches a MRI scan),musculoskeletal reshaping becomes more defined because one cannow use the pair of MRI images, namely the template and target,to aid with reshaping the template mesh [Gilles and Magnenat-Thalmann 2010; Gilles et al. 2006; Schmid et al. 2009]. The examplesin these papers demonstrate non-trivial musculoskeletal reshapinginvolving translation and spatially varying large rotations with alimited amount of volumetric stretching (Figure 14 in [Gilles andMagnenat-Thalmann 2010]). This is consistent with their choice ofthe similarity metric between the template and output shapes: theirreshaping energy tries to keep the distance of the output mesh tothe medial axis the same as the distance in the template [Gilles andMagnenat-Thalmann 2010], which biases the output against volumegrowth. For bones, a similar idea was also presented in [Schmidet al. 2011; Schmid and Magnenat-Thalmann 2008], where they didnot use a medial-axis term to establish similarity to a source mesh,but instead relied on a PCA prior on the shapes of bones, based ona database of 29 hip and femur bone shapes. Our work does notrequire any pre-existing database of shapes. Because our methoduses plasticity, it can accommodate large and spatially varying volu-metric stretching between the template and the subject. We do notneed a medical image for the template mesh. We only assume that , Vol. 1, No. 1, Article . Publication date: August 2020. • Bohan Wang, George Matcuk, and Jernej Barbič
Fig. 3.
Our shape deformation setup.
Our optimization discovers plasticstrains F p and vertex positions x so that the model is in elastic equilibriumunder the attachments, while meeting the medical image landmark andclosest point constraints as closely as possible. The presence of attachmentsis optional; our work handles both attached and un-attached objects. the template mesh is plausible. Of course, the template mesh itselfmight have been derived from or inspired by some MRI or CT scan,but there is no requirement that it matches any such scan. Anatomy Transfer.
Recently, great progress has been made on anatom-ically based simulations of humans. Anatomy transfer has been pio-neered by Dicko and colleagues [Dicko et al. 2013]. Anatomical mus-cle growth and shrinkage have been demonstrated in the “computa-tional bodybuilding” work [Saito et al. 2015]. Kadlecek et al. [2016]demonstrated how to transfer simulation-ready anatomy to a novelhuman, and Ichim et al. [Ichim et al. 2017] gave a state-of-the-artpipeline for anatomical simulation of human faces. Anatomy trans-fer and a modeling method such as ours are complementary becausethe former provides the means to interpolate known anatomies tonew subjects, whereas the latter provides a means to create theanatomies in the first place. Namely, anatomy transfer requires aquality anatomy template to serve as the source of anatomy transfer,which brings up the question of how one obtains such a template.Human anatomy is both extremely complex for each specific subject,and exhibits large variability in geometry across the population. Ac-curate templates can therefore only be created by matching them tomedical images. Even if one creates such a template, new templateswill always be needed to model the anatomical variability across theentire population; and this requires an anatomy modeling methodsuch as ours.
Given a template tet mesh of a soft tissue organ for a generic indi-vidual, as well as known optional attachments of the organ to otherobjects, our goal is to deform the tet mesh to match a medical imageof the organ of a new individual. We use the term “medical image”everywhere in this paper because this is standard terminology; thisdoes not refer to an actual 2D image, but to the 3D medical volume.We now describe how we mathematically model the attachmentsand medical image constraints.
We start with a template organ tet mesh M . We denote its vertexpositions by ¯x ∈ R n , where n is the number of tet mesh vertices.In this paper, we use bold text to represent the quantities for theentire mesh and non-bold text to represent the quantities for asingle vertex or element in the FEM mesh. We would like to dis-cover vertex positions x ∈ R n such that the organ shape obeysthe attachments to other organs (if they exist), and of course themedical image. The attachments are modeled by known materialpositions X i ∈ M , i = , . . . , t that have to be positioned at knownworld-coordinate positions x i ∈ R , i = , . . . t (Figure 3). Themedical image constraints come in two flavors. First, there are land-mark constraints whereby a point on a template organ is manuallycorresponded to a point in the medical image, based on anatomyknowledge. Namely, landmark constraints are modeled as materialpositions Y i ∈ M , i = , . . . , q , that are located at known world-coordinate positions y i ∈ R in the medical image. Observe thatlandmark constraints are mathematically similar to attachments.However, they have a different physical origin: attachments are aphysical constraint that is pulling the real-world organ to a knownlocation on another (fixed, un-optimized) object; for example, a mus-cle is attached to a bone. With landmarks, there is no such physicalforce in the real-world; namely, landmarks (and also closest-pointconstraints) are just medical image observations.The second type of medical image constraints are closest-point con-straints (“ICP markers”). They are given by known world-coordinatepositions z i ∈ R , i = , . . . r that have to lie on the surface of thedeformed tet mesh. Locations z i are easier to select in the medicalimage than the landmarks because there is no need to give any corre-spondence. As such, they require little or no medical knowledge, andcan be easily selected in large numbers. We simply went throughthe medical image slices and selected clear representative pointson the organ boundary. We then visually compared the templateand the target shape inferred by the ICP marker cloud. This guidedour positioning of the landmarks, which we place on anatomically“equal” positions in the template and the medical image. We con-sulted a medical doctor to help us interpret medical images, such asidentifying muscles in the scan, clarify ambiguous muscle bound-aries, placing difficult markers and attachments, and disambiguatingtendons. We model shape deformation using plastic deformation gradients,combined with a (small) amount of elastic deformation. In solidmechanics, plasticity is the tool to model large shape variations ofobjects, making it very suitable to model our desired shape defor-mation with large strains. Unlike using the elastic energy directly(without plasticity), plastic deformations have the advantage thatthey can arbitrarily and spatially non-uniformly and anisotropicallyscale the object. There is also no mathematical requirement that theyneed to respect volume preservation constraints. This makes plasticdeformations a powerful tool to model shapes. Our key idea is tofind a plastic deformation gradient F p at each tet of M , such thatthe FEM equilibrium shape under F p and any attachments matches , Vol. 1, No. 1, Article . Publication date: August 2020. odeling of Personalized Anatomy using Plastic Strains • 5 the medical image observations. Figure 3 illustrates our shape de-formation setting. In order to do so, we need to discuss the elasticenergy and forces in the presence of plastic deformations, whichwe do next.Plastic strain is given by a 3 × F p at each tetrahedronof M . For each specific deformed shape x ∈ R n , one can defineand compute the deformation gradient F between x and x at eachtet [Müller and Gross 2004]. The elastic deformation gradient F e can then be defined as F e = F F − p [Bargteil et al. 2007] (see Figure 5).Observe that for any shape x , there exists a corresponding plasticdeformation gradient F p such that x is the elastic equilibrium under F p ; namely F p = F . This means that the space of all plastic deforma-tion gradients F p is expressive enough to capture all shapes x . Theelastic energy of a single tet is defined as E ( F p , x ) = V ( F p ) ψ ( x , F p ) = V ( F p ) ψ (cid:0) F ( x ) F − p (cid:1) , (1)where V is the rest volume of the tet under the plastic deformation F p , and ψ is the elastic energy density function. We have V = | F p | V , where V is the tet’s volume in M , and | F p | is the determinant ofthe matrix F p . Elastic forces equal f e ( F p , x ) = d E ( F p , x )/ d x . Whensolving our optimization problem to compute F p in Section 3.4, wewill need the first and second derivatives of E ( F p , x ) with respectto x and F p . We provide a complete derivation of these terms inAppendix B.Our method supports any isotropic hyperelastic energy den-sity function ψ . In our examples, we use the isotropic stable neo-Hookean elastic energy [Smith et al. 2018], because we found it tobe stable and sufficient for our examples. Note that we do modelanisotropic plastic strains (and this is crucial for our method), sothat our models can stretch by different amounts in different di-rections. Observe that plastic strains are only determined up to arotation. Namely, let F p be a plastic strain (we assume det ( F p ) > F p = QS be the polar decompositionwhere Q is a rotation and S a 3 × symmetric matrix. Then, F p and S are the “same” plastic strain: the resulting elastic deformationgradients differ only by a rotation, and hence, due to isotropy of ψ , produce the same elastic energy and elastic forces. Note that itis not required that rotations Q match in any way at adjacent tets.We do not need to even guarantee that F p globally correspond toany specific “rest shape”, i.e., the F p are independent of each otherand may be inconsistent. This gives plastic deformation gradientmodeling a lot of flexibility. Hence, it is sufficient to model plasticstrains as symmetric 3 × F p as asymmetric matrix and parameterize it using a vector s ∈ R , F p = s s s s s s s s s . (2)We model plasticity globally using a vector s ∈ R m , where m is thenumber of tets in M . We note that Ichim et al. [2017] used such a6-dimensional parameterization to model facial muscle activations.In our work, we use it for general large-strain shape modeling. Ourapplication and optimization energies are different, e.g, Ichim etal. [2017] causes muscle shapes to follow a prescribed muscle firingfield, and biases principal stretches to be close to 1. Furthermore,we address the singularities arising with un-attached objects.
We now formulate our shape deformation problem. We first do sofor attached objects. An object is “attached” if there are sufficientattachment forces to remove all six rigid degrees of freedom, whichis generally satisfied if there are at least three attached non-colinearvertices. We find the organ’s shape that matches the attachment andthe medical image constraints by finding a plastic strain F p at eachtet, as well as static equilibrium tet mesh vertex positions x underthe attachments and plastic strain F p , so that the medical imageobservations are met as closely as possible,arg min s , x || L s || + α E MI ( x ) + β E a ( x ) , (3)subject to: f e ( F p ( s ) , x ) + f a ( x ) = , (4)where α ≥ β ≥ L is the plastic strain Laplacian . We define L as essentially the tet-meshLaplacian operator on the tets, 6-expanded to be able to operate onentries of s at each tet (precise definition is in Appendix A). TheLaplacian term enforces the smoothness of F p , i.e., F p in adjacenttets should be similar to each other. The second equation enforcesthe elastic equilibrium of the model under plastic strains F p and un-der the attachment forces f a . As such, our output shapes are alwaysin static equilibrium under the plastic strains F p , and both this equi-librium shape x and F p are optimized together; this is the key aspectof our work. The first equation contains the smoothness and themedical image (MI) observations; we discuss the attachment energy E a in the next paragraph. The medical image energy measures howclosely x matches the medical image constraints, E MI ( x ) = q (cid:213) i = || S x − y i || + r (cid:213) i = || z i − closestPoint ( x , z i )|| , (5)where S is the interpolation matrix that selects Y i , namely S ¯ x = Y . The function closestPoint ( x , z i ) ∈ R computes the closest point to z i ∈ R on the surface of the tet mesh with vertex positions x . Our treatment of attachments in Equations 3 and 4 deserves a spe-cial notice. Equation 4 is consistent with our setup: we are trying toexplain the medical images by saying that the organ has undergonea plastic deformation due to the variation between the template andcaptured individual. The shape observed in the medical image isdue to this plastic deformation and the attachments. We formulateattachment forces in Equation 4 as a “soft” constraint, i.e., f a ( x ) is modeled as (relatively stiff) springs pulling the attached organpoints to their position on the external object. This soft constraintcould in principle be replaced for a hard constraint where the at-tached positions are enforced exactly. We use soft constraints inour examples because they provide additional control to balanceattachments against medical image landmarks and ICP markers.These inputs are always somewhat inconsistent because it is impos-sible to place them at perfectly correct anatomical locations, dueto medical imaging errors. Hence, it is useful to have some leewayin adjusting the trade-off between satisfying each constraint type.With soft constraints, it is important to keep the spring coefficientin f a ( x ) high so that constraints are met very closely (under 0 . E a , we initially tried solving theoptimization problem of Equations 3 and 4 without it. This seems , Vol. 1, No. 1, Article . Publication date: August 2020. • Bohan Wang, George Matcuk, and Jernej Barbič Fig. 4.
Comparison to augmented deformation transfer.
The beam’s attachments (red) cause the beam to bend, whereas the ICP markers (blue) cause itto stretch 2x in one of the two transverse directions. Our method can easily recover such a shape deformation, whereas deformation transfer [Sumner andPopović 2004] cannot, even if augmented with an elastic energy.Fig. 5.
Plastic and elastic deformation gradient for a single tet. natural, but actually did not work. Namely, without E a , there isnothing in Equations 3 and 4 that forces the plastic strains to rea-sonable values. The optimizer is free to set F p to arbitrarily extremevalues, and then find a static equilibrium x under the attachmentforces. In our outputs, we would see smooth nearly tet-collapsingplastic strains that result in a static equilibrium x whereby the med-ical image constraints were nearly perfectly satisfied. Obviously,this is not a desired outcome. Our first idea was to add a term thatpenalizes the elastic energy E ( F p , x ) to Equation 3. Although thisworked in simple cases, it makes the expression in Equation 3 gen-erally nonlinear. Instead, we opted for a simpler and more easilycomputable alternative, namely add the elastic spring energy of allattachments, E a . This keeps the expression in Equation 3 quadraticin x and F p , which we exploit in Section 3.4 for speed. Observe that E a behaves similarly to the elastic energy: if the plastic strain causesa rest shape that is far from the attachment targets, then both E a and the elastic energy will need to “work” to bring the shape x to itstarget attachments. Similarly, if the plastic strain already did mostof the work and brought the organ close to its target, then neither E a nor the elastic energy will need to activate much.Because our units are meters and we aim to satisfy constraintsclosely, we typically use weights close to α = and β = inour examples. The weights α and β permit adjusting the trade-offbetween three desiderata: make plastic strains smooth, meet medicalimage observations, and avoid using too much elastic energy (i.e.,prefer to resolve shapes with plastic strains).Finally, we note that our formulation is different to approachesthat optimize the deformation gradient F directly (i.e., without anintermediary quantity such as the plastic deformation gradient). InFigures 4 and 7, we compare to two such approaches: deformationtransfer [Sumner and Popović 2004] and variational shape model-ing [Botsch and Kobbelt 2004]. We demonstrate that our methodbetter captures shapes defined using our inputs (landmarks, ICP markers, large spatially varying strains). Among all compared ap-proaches, the variational method in Figure 7 came closest to meetingour constraints, but there is still a visual difference to our method.We provide a further comparison to variational methods in Sec-tion 4.5. We adapt the Gauss-Newton method [Sifakis et al. 2005] to efficientlysolve the optimization problem of Equations 3 and 4 (example out-put shown in Figure 6). Doing so is not straightforward becausea direct application of the Gauss-Newton method results in largedense matrices that are costly to compute and store, causing themethod to fail on complex examples. Below, we demonstrate howto avoid these issues, producing a robust method capable of han-dling complex spatially varying plastic strains. Before settling onour specific Gauss-Newton approach, we also attempted to use theinterior-point optimizer available in the state-of-the-art Knitro op-timization library [Artelys 2019]. This did not work well becauseour problem is highly nonlinear. The interior-point method (IPM)worked well on simple examples, but was slow and not convergenton complex examples. IPM fails because it requires the constraintHessian, which is not available. When we approximated it, IPM gen-erated intermediate states too far from the constraints, and failed.The strength of our Gauss-Newton approach is that we only needconstraint gradients. Our method inherits the convergence proper-ties of the Gauss-Newton method. While not guaranteed to be locallyconvergent, Gauss-Newton is widely used because its convergencecan approach quadratic when close to the solution.We note that our method is designed for sparse medical imagelandmarks and ICP markers. In Figure 8, we give a comparison to arelated method from medical imaging which used an elastic energy,but with dense correspondences. Our method can produce a qualityshape even under sparse inputs, and can consequently work evenwith coarser MRI scans (such as our hip bone example; Figure 13).The ability to work with sparse markers also translates to lowermanual processing time to select the markers in the medical image.Because the object is attached, Equation 4 implicitly defines x as a function of s . The Gauss-Newton method uses the Jacobian J = d x / d s , which models the change in the static equilibrium x asone changes the plastic deformation gradient. It eventually relies on , Vol. 1, No. 1, Article . Publication date: August 2020. odeling of Personalized Anatomy using Plastic Strains • 7 Fig. 6.
17 muscles of the human hand extracted from MRI.
Observe that the template hand is bigger than the scanned hand. Pose is also different. Ourmethod solves this using bone attachments.Fig. 7.
Comparison to variational shape modeling.
In a variationalmethod [Botsch and Kobbelt 2004], the wiggles increase if one imposesa stricter constraint satisfaction. First row: under a small number of land-marks, variational methods with k = , produce a smooth and reasonableresult, albeit somewhat smoothening the rest shape. Middle row: undermore landmarks, it becomes more difficult for variational methods to meetthe landmarks while avoiding the wiggles. Bottom row: variational methodsproduce wavy results. Our method meets the landmarks and produces fewerwiggles. This is because the plastic deformation field can arbitrarily rotateand non-uniformly scale to adapt to the inputs; the elastic energy thenfinally irons out the kinks. the derivative of elastic forces with respect to the plastic deformationgradient, which we give in Appendix B. Because E MI ( x ) and E a ( x ) are quadratic functions of x , we rewrite Equations 3 and 4 asarg min x , s || Ls || + q + r + t (cid:213) k = c k || A k x + b k || , (6)s.t. f net ( s , x ) = , (7)where f net ( s , x ) = f e (cid:0) F p ( s ) , x (cid:1) + f a ( x ) is the net force on the mesh,and constant matrices, vectors and scalars A k , b k , c k are indepen-dent of s and x (we give them in Appendix F). The integer t denotesthe number of attachments. We now re-write Equations 6 and 7 sothat the plastic strains are expressed as s + ∆ s , and the equilibrium x as x + ∆ x , where ∆ x = J ∆ s . At iteration i of our Gauss-Newtonmethod, given the previous iterates x i and s i , we minimize a non-linearly constrained problem,arg min x i + , ∆ s i || L ( s i + ∆ s i )|| + q + r + t (cid:213) k = c k || A k (cid:0) x i + J ∆ s i (cid:1) + b k || , (8)s.t. f net ( s i + ∆ s i , x i + ) = . (9)After each iteration, we update s i + = s i + ∆ s i . Observe that Equa-tion 8 does not depend on x i + , and that the constraint of Equation 9is already differentially “baked” into Equation 8 via ∆ x = J ∆ s . Wetherefore first minimize Equation 8 for ∆ s i , using unconstrainedminimization; call the solution ∆ s i . A naive minimization requiressolving a large dense linear system of equations, which we avoid us-ing the technique presented at the end of this section. We regularize ∆ s i so that the corresponding F p is always positive-definite for eachtet; we do this by performing eigen-decomposition of the symmetricmatrix F p at each tet, and clamping any negative eigenvalues toa small positive value (we use 0.01). Our method typically did notneed to perform clamping in practice, and in fact such clamping isusually a sign that the method is numerically diverging, and shouldbe restarted with better parameter values.We then minimize the optimization problem of Equations 8 and 9using a 1D line search, using the search direction ∆ s i . Specifically,for η ≥ , we first solve Equation 7 with s ( η ) : = s i + η ∆ s i for x = x ( η ) using the Knitro library [Artelys 2019]. Direct solutionsusing a Newton-Raphson solver also worked, but we found Knitroto be faster. We then evaluate the objective of Equation 6 at x = x ( η ) and s = s ( η ) . We perform the 1D line search for the optimal η using , Vol. 1, No. 1, Article . Publication date: August 2020. • Bohan Wang, George Matcuk, and Jernej Barbič Fig. 8.
Elastic energy methods only work with dense markers.
In this figure, we compare to a state-of-the art medical imaging technique [Niculescuet al. 2009], whereby the output shape is calculated by minimizing an elastic energy of a template shape, subject to dense medical image markers. With densemarkers, elastic energy methods work well (left). As the constraints sparsify, elastic energy produces artifacts (middle). Our plasticity method (right) producesa good shape even with sparse markers. the Brent’s method [Press et al. 2007] because it does not require agradient.
Initial guess:
We solve our optimization problem by first assuminga constant s at each tet, starting from the template mesh as theinitial guess. This roughly positions, rotates and globally scales thetemplate mesh to match the medical image. We use the output asthe initial guess for our full optimization as described above. Fig. 9.
Convergence plots.
X-axis are iterations, and Y-axis is the optimiza-tion energy. The initial optimization energies are normalized to 1.0.
Optimization stages and stopping criteria:
We first do the opti-mization with attachments only. Upon convergence, we add thelandmarks, ignoring any ICP markers. This is because initially, themesh is far away from the target and the ICP closest locations areunreliable. After convergence, we disable the landmarks and enablethe ICP markers and continue optimizing. After this optimizationmeets a stopping criterium, we are done. Our output is thereforecomputed with ICP markers only; landmarks only serve to guidethe optimizer. This is because landmarks require a correct corre-spondence, and it is harder to mark this correspondence reliably in the scan than to simply select an ICP marker on the boundary ofan organ. We recompute the closest locations to ICP markers aftereach Gauss-Newton iteration. We stop the optimization if either ofthe following three criteria is satisfied: (i) reached the user-specifiedmaximal number of iterations (typically 20; but was as high as 80 inthe liver example), (ii) maximum error at ICP markers is less thana user-specified value (1mm for hand muscles), (iii) the progressin each iteration is too small, determined by checking if η is undera user-defined threshold (we use 0.01). Figure 9 shows the conver-gence of our optimization. Avoiding dense large linear systems.
Because the object is attached, ∂ f net ∂ x is square and invertible. Therefore, one can obtain a formulafor J by differentiating Equation 7 with respect to s , J = − (cid:16) ∂ f net ∂ x (cid:17) − ∂ f net ∂ s . (10)The matrix J is dense (dimensions 3 n × m ). Observe that becauseEquation 8 is quadratic in ∆ s i , minimizing it as done above to de-termine the search direction ∆ s i is equivalent to solving a linearsystem with the system matrix H , where H is the second derivative(Hessian matrix; dimension 6 m × m ) of Equation 8 with respectto ∆ s i . Because J is dense, H is likewise a dense matrix, H = L + J T (cid:0)(cid:213) k c k A k T A k (cid:1) J = L + Z T Z , (11)where Z = c A c A ... c r A r J = − c A c A ... c r A q + r + t (cid:16) ∂ f net ∂ x (cid:17) − ∂ f net ∂ s . (12)Therefore, when the number of tet mesh elements m is large, it isnot practically possible to compute H , store it explicitly in mem-ory or solve linear systems with it. To avoid this problem, we firsttried solving the system of equations using the Conjugate Gradient(CG) method. This worked, but was very slow (Table 1). The matrix Z ∈ R ( q + r + t ) × m is dense. In our complex examples, the numberof medical image constraints q + r + t is small (typically 10 - 800)compared to the dimension of s (6 m ; typically ~200,000). Our ideais to efficiently compute the solution to a system H y = h for any , Vol. 1, No. 1, Article . Publication date: August 2020. odeling of Personalized Anatomy using Plastic Strains • 9 Table 1.
Solving a single linear system of equations with H, using theconjugate gradients and our method. The naive direct solver failed in allcases. Note that H is a dense m × m matrix. The column t prep gives acommon pre-processing time for both CG and our method. Example 6 m ( q + r + t ) t prep CG OursHand muscle 237,954 1,143 17.5s 897.5s 9.5sHip bone 172,440 1,497 11.4s 408.9s 7.2sLiver 259,326 1,272 24.4s 1486.7s 10.0sright-hand side h using the Woodbury matrix identity [Woodbury1950], where we view L as a “base” matrix and Z T Z a low-rankperturbation. Before we can apply Woodbury’s identity, we needto ensure that the base matrix is invertible. As we prove in Appen-dix A, the plastic strain Laplacian L is singular with six orthonormalvectors ψ i in its nullspace (assuming that M is connected). Each ψ i is a vector of all ones in component i of s , i = , . . . , √ m for normalization. It follows fromthe Singular Lemma (i) (Section 3.5) that L is also singular with thesame nullspace vectors. Therefore, we decompose H = (cid:16) L − (cid:213) i = ψ i ψ i T (cid:17) + (cid:16) Z T Z + (cid:213) i = ψ i ψ i T (cid:17) = B + ˆZ T ˆZ , (13)where B = L − (cid:205) i = ψ i ψ i T and ˆZ is matrix Z with an additionaladded 6 added rows ψ i T . By the Singular Lemma (iii) (Section 3.5), B is now invertible, and we can use Woodbury’s identity to solve y = H − h = (cid:16) B − − B − ˆZ T (cid:16) I + ˆZB − ˆZ T (cid:17) − ˆZB − (cid:17) h . (14)We rapidly compute Z , without ever computing or forming J , bysolving sparse systems ∂ f net ∂ x z = c k A k T , for k = , . . . , q + r + t . Observe that this sparse system matrix is symmetric and the samefor all k . We factor it once using the Pardiso solver and then solve themultiple right-hand sides in parallel. The matrix B is constant, andwe only need to factor it once for the entire optimization. Finally, thematrix I + ˆZB − ˆZ T ∈ R ( q + r + t ) × ( q + r + t ) is small, and so invertingit is fast. We analyze the performance of our algorithm in Table 1. Fig. 10.
Illustration of theSingular Lemma.
In this paper, there are two occasionswhere we have to solve a singular sparselinear system with known nullspacevectors. Such systems occur often inmodeling of un-attached objects, e.g.,finding static equilibria, solving Laplaceequations on the object’s mesh, an-imating with rotation-strain coordi-nates [Huang et al. 2011], or comput-ing modal derivatives [Barbič and James2005]. Previous work solved such systems ad-hoc, and the under-lying theory has not been stated or developed in any great de-tail. We hereby state and prove a lemma that comprehensivelysurveys the common situations arising with singular systems incomputer animation and simulation, and back the lemma with amathematical proof (Appendix C). Recall that the nullspace of amatrix A ∈ R p × p is N( A ) = { x ∈ R p ; Ax = } , and the range of A is R( A ) = { Ax ; x ∈ R p } . Both are linear vector subspaces of R p . Singular Lemma:
Let the square symmetric matrix A ∈ R p × p besingular with a known nullspace spanned by k linearly independentvectors ψ , . . . , ψ k . Then the following statements hold:(i) N( A ) and R( A ) are orthogonal. Every vector x ∈ R p can beuniquely expressed as x = n + r , where n ∈ N( A ) and r ∈ R( A ) . Vector r is orthogonal to n and to ψ i for all i = , . . . , k (Figure 10).(ii) Let b ∈ R( A ) . Then, the singular system Ax = b has a uniquesolution x with the property that x is orthogonal to ψ i for all i = , . . . , k . This solution can be found by solving the non-singular linear system
A ψ . . . ψ k ψ T . . . ... ... ... ψ Tk . . . xλ ... λ k = b ... . (15)All other solutions equal x + (cid:205) ki = µ i ψ i for some scalars µ i ∈ R . (iii) For any scalars α i (cid:44) , the matrix B = A + (cid:205) ki = α i ψ i ψ Ti isinvertible. If ψ i are orthonormal vectors, then the solution to By = h equals y = x + (cid:205) ki = λ i α i ψ i , where x and λ i are solutions to Equa-tion 15 with b = proj R( A ) h = h − (cid:205) ki = ( ψ Ti h ) ψ i . We give the proofof the singular lemma in Appendix C.
The difficulty with un-attached objects is that we now have f net = f e , and the equation f e (cid:0) F p ( s ) , x (cid:1) = x for a fixed plastic state s . This can be intuitively easily understood:one can arbitrarily translate and rotate any elastic equilibrium shape x under the given plastic state s ; doing so produces another elasticequilibrium shape. The space of solutions x is 6-dimensional. Thismeans that we can no longer uniquely solve Equation 7 for x duringour line search of Section 3.4. Furthermore, the square tangentstiffness matrix K ( F p ( s ) , x ) = ∂ f e (cid:0) F p ( s ) , x (cid:1) ∂ x (16)is no longer full rank. In order to address this, we now state andprove the following Nullspace Lemma. Nullspace Lemma:
The nullspace of the tangent stiffness matrixof an elastoplastic deformable object in static equilibrium x underplasticity, is 6-dimensional. The six nullspace vectors are ψ i : = [ e i , e i , . . . , e i ] , where e i ∈ R is the i -th standard basis vector, and ψ + i : = [ e i × x , e i × x , . . . , e i × x n ] , for i = , , . To the best of our knowledge, this fact of elasto-plasto-staticshas not been stated or proven in prior work. It is very useful whenmodeling large-deformation elastoplasticity, as real objects are oftenun-attached, or attachments cannot be easily modeled. We give aproof in Appendix D. To accommodate un-attached objects, it istherefore necessary to stabilize the translation and rotation. Fortranslations, this could be achieved easily by fixing the position ofany chosen vertex. Matters are not so easy for rotations, however.Our idea is to constrain the centroid of all tet mesh vertices to aspecific given position t , and to constrain the “average rotation” of , Vol. 1, No. 1, Article . Publication date: August 2020. Fig. 11.
Un-attached optimization of a hip bone shape to a CT scan.
The scanned bone is smaller and has a substantially different shape to the template.Fig. 12.
Un-attached optimization of a liver shape to a CT scan.
Our method successfully captures the large shape variation between the template andthe scan. This figure also demonstrates that our method makes it possible to transfer the rendering textures and uv coordinates from the template onto theoutput.Fig. 13.
Attached optimization of a hip muscle (gluteus medius) to a MRI scan.
Our method successfully captures the large shape variation betweenthe template and the scan. the model to a specific given rotation R . We achieve this using thefamiliar “shape-matching” [Müller et al. 2005], by imposing that therotation in the polar decomposition of the global covariance matrixmust be R . We therefore solve the following optimization problem,arg min x , s , R , t || L s || + α E MI ( x ) (17)s.t. f e (cid:0) F p ( s ) , x (cid:1) = , (18) (cid:169)(cid:173)(cid:171)(cid:213) j ∈ D w j x j (cid:170)(cid:174)(cid:172) − (cid:213) j ∈ D w j X j = t , (19)Polar (cid:169)(cid:173)(cid:171)(cid:213) j ∈ D w j (cid:16) x j − t (cid:17) (cid:16) X j − (cid:0) (cid:213) k ∈ D w k X k (cid:1)(cid:17) T (cid:170)(cid:174)(cid:172) = R , (20)where D is the set of points on the mesh surface where we have ei-ther a landmark or an ICP constraint, w j is the weight of a point, X j is the position of vertex j in M and Polar ( F ) is the polar decomposi-tion function that extracts the rotational part of a matrix F . We set all weights equal, i.e., w j = /| D | . We choose the set D as opposedto all mesh vertices so that we can easily perform optimization withrespect to R and t (next paragraph). We assume that our argumentmatrices F to Polar are not inversions, i.e., det ( F ) > , which es-tablishes that Polar ( F ) is always a rotation and not a mirror. Thisrequirement was easily satisfied in our examples, and is essentiallydetermined by the medical imaging constraints; the case det ( F ) < x , s whilekeeping R , t fixed and vice-versa (Figure 11, 12). Rigid transforma-tions do not affect smoothness of s so we do not need to considerit when optimizing R , t . We need to perform two modifications toour Gauss-Newton iteration of Section 3.4. The first modificationis that we need to simultaneously solve Equations 18, 19 and 20when determining the static equilibrium in the current plastic state s . As with attached objects, we do this using the Knitro optimizer. In , Vol. 1, No. 1, Article . Publication date: August 2020. odeling of Personalized Anatomy using Plastic Strains • 11 order to do this, we need to compute the first and second derivativesof the stabilization constraints in Equations 19 and 20 (Section 3.7).The second modification is needed because the tangent stiffness ma-trix K ( F p ( s ) , x ) , as explained above, is now singular with a known6-dimensional nullspace. In order to compute the Jacobian matrix J using Equation 10, we use our Singular Lemma (ii) (Section 3.5).Note that the right-hand side is automatically in the range of K be-cause Equation 10 was obtained by differentiating a valid equation,hence Equation 10 must also be consistent. Polar ( F ) Previous work computed first and second-order time derivatives ofthe rotation matrix in polar decomposition [Barbič and Zhao 2011],or first derivative with respect to each individual entry of F [Chaoet al. 2010; Twigg and Kačić-Alesić 2010]. In our work, we need thefirst and second derivatives of R with respect to each individualentry of F . We found an elegant approach to compute them usingSylvester’s equation, as follows. Observe that Polar ( F ) = FS − , where S is the symmetric matrix in the polar decomposition. Becausedet ( F ) > , S is positive-definite and uniquely defined as S = √ F T F . To compute the first-order derivatives, we start from F = RS , anddifferentiate, ∂ F ∂ F i = ∂ R ∂ F i S + R ∂ S ∂ F i , hence ∂ R ∂ F i = (cid:18) ∂ F ∂ F i − R ∂ S ∂ F i (cid:19) S − , (21)Therefore, we need to compute ∂ S / ∂ F i . We have F T F = S , and thus ∂ F T F ∂ F i = ∂ S ∂ F i S + S ∂ S ∂ F i , (22)i.e., this is the classic Sylvester equation for the unknown matrix ∂ S ∂ F i [Sylvester 1884]. The Sylvester equation AX + X B = C can besolved as (cid:16) B T ⊕ A (cid:17) − vec ( X ) = vec ( C ) , (23)where ⊕ is the Kronecker sum of two matrices. In our case,vec ( ∂ S ∂ F i ) = ( S ⊕ S ) − vec ( ∂ F T F ∂ F i ) . (24)The computation of second-order derivatives follows the samerecipe: differentiate the polar decomposition and solve a Sylvesterequation. We give it in Appendix E.We can now compute the gradient and Hessian of our stabiliza-tion constraints. The translational constraint is linear in x and canbe expressed as W x − d = , where W is a 3 × n sparse matrix.Although Polar is not linear, the argument of Polar is linear in x . Therotational constraint can be expressed as Polar ( W x − d ) − ¯ R = , where W is a 9 × n sparse matrix. The Jacobian of the transla-tional constraint is W , and the Hessian is zero. For the rotationalconstraint, the Jacobian is ∂ R ∂ F : W ∈ R × n and the Hessian is ( W T : ∂ R ∂ F : W ) ∈ R × n × n , where : denotes tensor contraction. Our shape deformation setup is similar to standard shape deforma-tion problems in computer graphics. In fact, we first attempted tosolve the shape deformation problem with as-rigid-as-possible en-ergy (ARAP) [Sorkine and Alexa 2007], bounded biharmonic weights
Fig. 14.
Standard volume-based shape deformation methods resultin wiggly and spiky artefacts.
The shown hand Palmar interossei musclehas a tendon on one side and no tendon on the other; both ends are attachedto a bone (light blue). Definitions of the acronyms are in Section 3.8. TheMRI landmark constraints are shown in dark blue. The shape deformationbetween the template muscle and the shape in MRI consists of large andspatially varying stretches. Our method successfully models this deforma-tion. We note that spikes cannot be avoided by using, say, a spherical regionfor the constraints as opposed to a point; the non-smoothness just movesto the spherical region boundary.Fig. 15.
Comparison between a surface energy and volumetric en-ergy.
In both examples (bunny and hand), we performed non-rigid iterativeclosest point alignment between a template triangle mesh, and a collectionof target ICP markers (13 for bunny and 456 for the hand). For the bunny,we manually placed the markers to greatly enlarge bunny’s left ear. For thehand, we placed the markers on a pre-existing target hand surface meshthat has different geometric proportions and mesh topology as the template.Template is a man and target is a woman. We then solved the ICP problemusing the surface-based ARAP energy, volume-based ARAP energy andour volumetric plastic strains. Our method produces smooth artifact-freeoutputs. (BBW) [Jacobson et al. 2011], biharmonic weights with linear preci-sion (LBW) [Wang et al. 2015a], and a Finite Element Method staticsolver (FEM) [Barbič et al. 2009]. Unfortunately, none of the meth-ods worked well. Figures 14 and 15 demonstrate that these methodsproduce non-smooth shapes with spikes (ARAP, BBW, FEM), orwiggles (LBW).Mathematically, the reason for the spikes in ARAP, BBW and FEMis that point constraints for second-order methods are inherently , Vol. 1, No. 1, Article . Publication date: August 2020. flawed. As one refines the solution by adding more tetrahedra, thesolution approaches a spiky point function at each point constraint,which is obviously not desirable. This mathematical issue is exposedin our work because our shape deformation consists of large spatiallyvarying stretches. Often, the template mesh needs to be stretched ∼
2x or more along some (or several) coordinate axes. The medicalimage constraints are distributed all around the muscle, pulling indifferent directions and essentially requesting the object to undergoa spatially non-uniform and anisotropic scale. This exacerbates thespikiness for second-order methods. We note that these problemscannot be avoided simply by using an elastic energy that permitsvolume growth. Namely, Drucker’s stability condition [Drucker1957] requires a monotonic elastic energy increase with increase instrain. An elastic energy therefore must penalize strain increases ifit is to be stable; and this impedes large-strain modeling in methodsthat rely purely on an elastic energy. Our plasticity method doesnot penalize large strains and thus avoids this problem. Spikes canbe avoided by using a higher-order variational method such asLBW. However, our experiments indicate that such methods sufferfrom wiggles when applied to medical imaging problems (see alsoFigures 7 and 23).
We extracted muscles of the human hand and the hip muscle froman MRI scan, and a hip bone and a liver from a CT scan. We analyzethe performance of our method in Table 2. In Figure 16, we givehistograms of the magnitude of the difference between the positionsof the medical image markers and their output positions. It canbe seen that our method produces shapes that generally matchthe medical image constraints to 0.5mm or better. In Figure 17, wedemonstrate that the output quality of our tetrahedra is still good; ifneeded, this could be further improved by re-meshing [Bargteil et al.2007]. Figures 18 and 19 superimpose our output meshes on the CTand MRI scans, respectively. In Figure 20, we compare to a recentimplicit point set surface reconstruction method [Huang et al. 2019].In Figure 21, we evaluate our method in the presence of knownground truth plastic deformation gradients. Figure 22 provides acomparison to surface-based methods.
In our muscle hand example, we extracted 17 hand muscles from anMRI scan (Figure 6). We obtained the scan and the already extractedbone meshes from [Wang et al. 2019]; scan resolution is 0.5mm x0.5mm x 0.5mm . We considered two “templates”, the first one fromthe Centre for Anatomy and Human Identification at the Universityof Dundee, Scotland [2019], and the second one from Zygote [Zygote2016]. We used the first one (Figure 6, left) because we found it to bemore medically accurate (muscles insert to correct bones). Muscleanatomy of a human hand is challenging (Figure 6). We model allmuscle groups of the hand, namely the thenar eminence (thumb),hypothenar eminence (below little finger), interossei muscles (pal-mar and dorsal) (between metacarpal bones), adductor pollicis (softtissue next to the thumb, actuating thumb motion), and lumbricals(on the side of the fingers at the fingers base). Our template modelsthe correct number and general location of the muscles, but there
Fig. 16.
Output ICP error histograms.
Each medical image marker con-tributes 1 entry to the histogram. The hand muscles histogram is a combinedhistogram for all the 17 hand muscles.Fig. 17.
Minimum tetrahedral dihedral angles before and after ouroptimization.
It can be seen that the output angles are still relatively large.As expected, the output angles are somewhat worse than the initial ones,as the object has undergone a plastic deformation. are large muscle shape differences between the template subjectand the scanned subject (Figure 1). We solve the optimization prob-lem of Equations 3 and 4 separately for each muscle, starting fromthe template mesh as the initial guess. In our results, this producesmuscles that match the attachments and medical image constraintsmarkers at 0.5 mm or better, which is at, or better than, the accuracyof the MRI scanner.
During pre-processing, wemanually mark as many reliable points as possible on the boundaryof each muscle ( ∼ −
20 landmarks and ∼ −
100 ICP markersper muscle) in the MRI scans. This process took approximately 5minutes per muscle.
The template muscles are modeledas triangle meshes. We build a tetrahedral mesh for each muscle. , Vol. 1, No. 1, Article . Publication date: August 2020. odeling of Personalized Anatomy using Plastic Strains • 13
Table 2.
The statistics for our examples: M ; e init = error between our template mesh and theICP constraints; e final = error between our result and the ICP markers. The first and second reported error numbers are the average and maximum error,respectively. In the hand example, there are 17 groups of muscles; “min”, “med” and “max” refers to the smallest, representative median and largest musclegroup; “max-m” is the example with the largest number of ICP constraints. Example e init [ mm ] e final [ mm ] Hand muscle (min) 4,912 20,630 15 8 14.2 yes 0.53 / 2.22 0.06 / 0.14Hand muscle (med) 6,260 31,737 32 12 12.8 yes 0.62 / 1.56 0.11 / 0.55Hand muscle (max) 8,951 42,969 96 11 14.2 yes 3.35 / 12.82 0.11 / 0.34Hand muscle (max-m) 7,552 34,966 151 18 20.3 yes 3.28 / 9.11 0.16 / 0.47Hip muscle (Fig 13) 6,793 34,118 82 21 28.3 yes 7.41 / 21.27 0.39 / 1.85Hip bone 6,796 28,740 499 34 49.2 no 4.12 / 14.27 0.25 / 1.30Liver 11,392 43,221 424 80 128.3 no 9.00 / 33.68 0.21 / 4.81Hand surface (Fig 15) 11,829 49,751 456 31 43.8 no 4.87 / 16.78 0.07 / 0.86
Fig. 18.
Our extracted organs match the medical image.
The intersec-tion of the output mesh with the medical image slices is shown in green.
Our tet meshes conform to the muscle’s surface triangle mesh; thisrequirement could be relaxed. For each muscle in the template, weattach its tet mesh to the bones using soft constraints. We do thisby marking where on one or multiple bones this muscle inserts;to do so, we consulted a medical doctor with specific expertisein anatomy. For each bone triangle mesh vertex that participatesin the insertion, we determine the material coordinates (i.e., tetbarycentric coordinates) in the closest muscle tet. We then form asoft constraint whereby this muscle material point is linked to thebone vertex position using a spring.
We note that we have alsoattempted to model the muscle shapes directly using segmentation,simply from an MRI scan. Recent work has demonstrated that thiscan be done for hand bones [Wang et al. 2019], and we attempteda similar segmentation approach for muscles. However, given thatthe muscles touch each other in many places (unlike bones), thecontrast in the MRI scan was simply not sufficient to discern theindividual muscles. Our conclusion is that a segmentation approachis not feasible for hand muscles, and one must use a pre-existinganatomically accurate template as in our work.
Many hand musclesare in close proximity to one another and several are in continuouscontact. One strategy to resolve contact would be to incorporate con-tact into our optimization (Equations 3 and 4). This approach is notvery practical, as unilateral constraints are very difficult to optimize.Furthermore, such an approach couples all (or most) muscles, andrequires one to solve an optimization problem with a much largernumber of degrees of freedom. It is extremely slow at our resolu-tion; it overpowers our machine. Instead, we optimize each muscleseparately. Of course, the different muscles inter-penetrate eachother, which we resolve as follows. For each muscle, we already po-sitioned the MRI constraints so that they are at the muscle boundary.Therefore, observe that if our marker positioning and the solutionto the optimization problem were perfect, inter-penetration wouldbe minimal or non-existent. The marker placement is relativelystraightforward at the boundary between a muscle and anothertissue (bone, fat, etc.) due to good contrast. However, placing mark-ers at the boundary between two continuously colliding muscles , Vol. 1, No. 1, Article . Publication date: August 2020.
Fig. 19.
Interpenetration removal.
The yellow lines are muscle cross-sections in this representative MRI slice of the hand. It can be seen that ourinterpenetration removal method successfully removes penetrations, with-out modifying the interpenetration-free sections of the muscles’ boundary. is less precise, due to a lower MRI contrast between adjacent mus-cles. This is the main cause of the inter-penetrations. We removethe inter-penetrations with FEM simulation because it producessmooth organic shape changes; note that alternatively, geometricapproaches could also be used [Schmid et al. 2009]. Specifically inour work, for a pair of inter-penetrating muscles, we run collisiondetection to determine the set of triangles of each muscle that areinside the volume of the other muscle. On each muscle, we thendetermine the set of tetrahedra that are adjacent to the collisionarea. We then slightly enlarge this set, by including any tet within a5-ring of tets. We then run a FEM contact simulation just on thesetwo tetrahedral sets on the two muscles. The FEM simulation pushesthe muscle surface boundaries apart, without displacing the restof the muscle (Figure 19). We handle contact islands of multiplemuscles by running the above procedure on two muscles, then for athird muscle against the first two muscles, then the fourth againstthe first three, and so on.
In our second anatomical example (Figure 11), we apply our methodto a CT scan of the human right hip bone (pelvis). We obtained thetemplate from the human anatomy model of Ziva Dynamics [ZivaDynamics 2019], and the CT scan from the “KidneyFullBody” medi-cal image repository [Stephcavs 2019]. The template and the scannedhip bone differ substantially in shape, and this is successfully cap-tured by our method.
In our third anatomical example, we apply our method to a CTscan of the human liver (Figure 12). We purchased a textured livertriangle mesh on TurboSquid [Turbosquid 2019]. We subdivided itand created a tet mesh using TetGen [Hang Si 2011]. This serves as our “template”. We used a liver CT scan from the “CHAOS” medicalimage repository [Kavur et al. 2019]. We then executed our methodto reshape the template tet mesh to match the CT scan. Much likewith the hip bone, our method successfully models the large dif-ferences between the template and the scanned liver. Finally, weembedded the TurboSquid triangle mesh into the template tet mesh,and transformed it with the shape deformation of the tet mesh. Thisproduced a textured liver mesh (Figure 12) that matches the CTscan.
In our fourth anatomical example, we apply our method to a MRIscan of a female human hip muscle (gluteus medius) (Figure 13). Weobtained the data from The Cancer Imaging Archive (TCIA) [Clarket al. 2013]. The image resolution is 384 × ×
240 with voxel spac-ing of 1mm, which is 2x coarser to the hand MRI dataset. We usethe template mesh from the human anatomy model of Ziva Dynam-ics [Ziva Dynamics 2019]. We subdivided it and created a tet meshfor it using TetGen [Hang Si 2011]. Because the muscle is attachedto the hip bone and the leg bone, we needed to first extract the bonesfrom the MRI scan; we followed the method described in [Wanget al. 2019]. Note that the subject in the Ziva Dynamics template ismale. The template and the scanned hip muscle differ substantiallyin shape, and this is successfully captured by our method. (a) Our output (b) [Huang et. al 2019]
Fig. 20.
Comparison to [Huang et al. 2019].
Top row: hip bone. Bottomrow: liver. Red points are the markers from the CT scan. We used the publiclyavailable implementation of [Huang et al. 2019] to compute the normals,followed by screened Poisson surface reconstruction using the points andcomputed normals [Kazhdan and Hoppe 2013]. We used this combinationbecause it produced better results than running [Huang et al. 2019] directly.It can be seen that our method produces shapes that match the groundtruth data more closely.
We compare our method to variational shape modeling methods onan illustrative 1D example. Note that Figure 7 gave a comparisonon 3D muscle geometry. Consider an elastic 1D line segment whoseneutral shape is the interval [ , ] , and study its longitudinal 1Ddeformation under the following setup. Let us prescribe hard attach-ments whereby we attach endpoint 0 to position 0 , and endpoint 1 , Vol. 1, No. 1, Article . Publication date: August 2020. odeling of Personalized Anatomy using Plastic Strains • 15 Fig. 21.
Quantitative evaluation on ground truth.
We first performeda dynamic FEM simulation with plasticity [Irving et al. 2004], wherebythe back of the dragon is fixed (red), and the head was pulled to the leftwith uniform force. This produced our ground truth plastic deformationgradients. We then selected 488 sparse triangle mesh vertices as landmarksand ICP markers, and ran our method to compute F p and shape x . It canbe seen that F p and x match the ground truth closely.Fig. 22. Comparison to surface-based methods.
We compare to the“PRIMO” method [Botsch et al. 2006] because this method is a good rep-resentative choice. It also has fewest artifacts in Figure 10 of the shapedeformation survey [Botsch and Sorkine 2008]. Our methods produces aclearly superior result in the challenging scenario where the ICP markerswere placed to anisotropically stretch the beam in one transverse direction.Our method also passes the standard shape deformation benchmark [Botschand Sorkine 2008]. to position 2 . Furthermore, assume landmarks whereby point 1 / / , and point 1 / / . Effectively in this setup, we arespecifying that the subintervals [ , / ] and [ / , ] do not stretch,whereas the subinterval [ / , / ] stretches from its original lengthof 1 / / , (5x stretch). A variational formulation of order r is,min x ( t )∈C r ∫ (cid:0) d r xdt (cid:1) dt + α (cid:16)(cid:0) x ( ) − (cid:1) + (cid:0) x ( ) − (cid:1) (cid:17) , (25)where C r denotes all functions [ , ] → R whose derivatives ex-ist and are continuous up to order r . We solved these problemsanalytically in Mathematica for r = , , , each time for three rep-resentative values α , and compared them (see Figure 23) to our method in 1D,min f p , x ( t )∈C r ∫ (cid:220) f p dt + α (cid:16)(cid:0) x ( ) − (cid:1) + (cid:0) x ( ) − (cid:1) (cid:17) + β ∫ (cid:0) (cid:219) xf p − (cid:1) dt (26)s. t. x ( t ) = arg min x ( t )∈C r ∫ (cid:0) (cid:219) xf p − (cid:1) dt . (27)Observe that, in the same vein as in 3D, we can decompose (cid:219) x = f e f p , whereby scalar functions f e and f p and the β term are theequivalents of the elastic and plastic deformation gradients, andthe elastic energy, respectively. It can be seen that our methodproduces a better fit to the data and a significantly less wigglysolution, compared to variational methods (Figure 23). Fig. 23.
Comparison to variational methods.
A 1D string of length 1 isattached on both ends. The left attachment is fixed. The right attachmentis moved to coordinate 2. The spring thus stretches longitudinally whiletrying to obey the two landmarks at t = / and t = / . Y-axis gives thedeformed 1D position of the corresponding material point on the X-axis.Big and small dots denote the attachments and landmarks, respectively.For each method, we plot the result under a few representative parameters.With variational methods, one has to either give up meeting the landmarks,or the curve becomes wiggly. Similar behavior can also be observed in 3Dvariational results (Figure 7). Our method produces a deformation profilewhereby the slope (total 1D deformation gradient) on all three subintervals [ , / ] , [ / , / ] , [ / , ] approximately matches the slope implied by theattachments and landmarks (1, 5, 1, respectively). Note that the shownoutput curves are only C r − and not in C r at the two juncture points / , / however, their r -th derivative is integrable and they are the optimalCauchy limit of a score-decreasing sequence of curves in C r . We gave a shape deformation method that can model objects under-going large spatially varying strains. Our method works by comput-ing a plastic deformation gradient at each tet, such that the meshdeformed by these plastic deformation gradients matches the pro-vided sparse landmarks and closest-point constraints. We supportboth constrained and unconstrained objects; the latter are supportedby giving a numerical method to solve large sparse linear systemsof equations with a known nullspace. We applied our method to theextraction of shapes of organs from medical images. Our method has , Vol. 1, No. 1, Article . Publication date: August 2020. been designed to extract as much information as possible from MRIimages, despite the soft boundaries between the different organs.We extracted hand muscles, a liver, a hip bone and a hip muscle.Our method does not require a dense target mesh; only a sparseset of observations is needed. If a dense target mesh is available, theproblem becomes somewhat easier, as one can then use standardICP algorithms. However, medical images contain many ambigui-ties and regions where there is not a sufficient contrast to clearlydisambiguate two adjacent medical organs; making it impractical toextract dense target meshes. We apply our method to solid objects,but our plastic strain shape deformation method could also be usedfor shells (cloth). Doing so would require formulating the elasticenergy of a plastically deformed FEM cloth, and computing the en-ergy gradients and Hessians with respect to the plastic parameters.The size of the small square dense matrix that we need to invert inour incremental solve is three times the number of markers. Whilewe easily employed up to a thousand of markers in our work, ourmethod will slow down under a large number of markers. We donot re-mesh our tet meshes during the optimization. If the plasticstrain causes some tetrahedra to collapse or nearly collapse, thiswill introduce numerical instabilities. Although not a problem inour examples (see Figure 17), such situations could be handled by re-meshing the tet mesh during the optimization [Bargteil et al. 2007].Our method requires a plausible template with a non-degenerate tetmesh. Re-meshing is important future work as it could extend thereach of our method, enabling one to start the optimization simplyfrom a sphere tet mesh.
REFERENCES
M. Alexa, A. Angelidis, M.-P. Cani, S. Frisken, K. Singh, S. Schkolne, and D. Zorin. 2006.Interactive Shape Modeling. In
ACM SIGGRAPH 2006 Courses
ACM Trans. on Graphics
28, 3 (2009).J. Barbič and D. L. James. 2005. Real-time subspace integration for St. Venant-Kirchhoffdeformable models.
ACM Trans. on Graphics
24, 3 (2005), 982–990.J. Barbič and Y. Zhao. 2011. Real-time Large-deformation Substructuring.
ACM Trans.on Graphics (SIGGRAPH 2011)
30, 4 (2011), 91:1–91:7.A. W. Bargteil, C. Wojtan, J. K. Hodgins, and G. Turk. 2007. A finite element methodfor animating large viscoplastic flow. In
ACM Transactions on Graphics (SIGGRAPH2007) , Vol. 26. 16.Miklós Bergou, Saurabh Mathur, Max Wardetzky, and Eitan Grinspun. 2007. TRACKS:Toward directable thin shells.
ACM Trans. on Graphics (SIGGRAPH 2007)
26, 3 (2007),50:1–50:10.Mario Botsch and Leif Kobbelt. 2004. An Intuitive Framework for Real-Time FreeformModeling.
ACM Trans. on Graphics (SIGGRAPH 2004)
23, 3, 630–634.Mario Botsch, Mark Pauly, Markus Gross, and Leif Kobbelt. 2006. PriMo: CoupledPrisms for Intuitive Surface Modeling. In
Eurographics Symp. on Geometry Processing .11–20.M. Botsch and O. Sorkine. 2008. On linear variational surface deformation methods.
IEEE Trans. on Vis. and Computer Graphics
14, 1 (2008), 213–230.C. Erolin. 2019. Hand Anatomy. University of Dundee, Centre for Anatomy and HumanIdentification. https://sketchfab.com/anatomy_dundee/collections/hand-anatomy.I. Chao, U. Pinkall, P. Sanan, and P. Schröder. 2010. A Simple Geometric Model forElastic Deformations.
ACM Transactions on Graphics
29, 3 (2010), 38:1–38:6.W. Chen, F. Zhu, J. Zhao, S. Li, and G. Wang. 2018. Peridynamics-Based FractureAnimation for Elastoplastic Solids. In
Computer Graphics Forum , Vol. 37. 112–124.K. Clark, B. Vendt, K. Smith, J. Freymann, J. Kirby, P. Koppel, S. Moore, S. Phillips,D. Maffitt, M. Pringle, L. Tarbox, and F. Prior. 2013. The Cancer Imaging Archive(TCIA): maintaining and operating a public information repository.
J Digit Imaging
26, 6 (Dec 2013), 1045–1057.Ali Hamadi Dicko, Tiantian Liu, Benjamin Gilles, Ladislav Kavan, Francois Faure, OlivierPalombi, and Marie-Paule Cani. 2013. Anatomy Transfer.
ACM Trans. on Graphics(SIGGRAPH 2013)
32, 6 (2013), 188:1–188:8.Daniel Charles Drucker. 1957.
A definition of stable inelastic material . Technical Report.DTIC Document. B. Gilles and N. Magnenat-Thalmann. 2010. Musculoskeletal MRI segmentation usingmulti-resolution simplex meshes with medial representations.
Med. Image Anal.
Medical Image Computingand Computer-Assisted Intervention – MICCAI 2006 . 289–296.Benjamin Gilles, Lionel Reveret, and Dinesh Pai. 2010. Creating and animating subject-specific anatomical models.
Computer Graphics Forum
29, 8 (2010), 2340–2351.Leo Grady. 2006. Random walks for image segmentation.
IEEE Trans. on Pattern Analysisand Machine Intelligence
28, 11 (2006), 1768–1783.Hang Si. 2011. TetGen: A Quality Tetrahedral Mesh Generator and a 3D DelaunayTriangulator.Jin Huang, Yiying Tong, Kun Zhou, Hujun Bao, and Mathieu Desbrun. 2011. InteractiveShape Interpolation through Controllable Dynamic Deformation.
IEEE Trans. onVisualization and Computer Graphics
17, 7 (2011), 983–992.Z. Huang, N. A. Carr, and T. Ju. 2019. Variational implicit point set surfaces.
ACMTrans. on Graphics (SIGGRAPH 2019)
38, 4 (2019).A.E. Ichim, P. Kadlecek, L. Kavan, and M. Pauly. 2017. Phace: Physics-based FaceModeling and Animation.
ACM Trans. on Graphics (SIGGRAPH 2017)
36, 4 (2017).T. Igarashi, T. Moscovich, and J. F. Hughes. 2005. As-rigid-as-possible shape manipula-tion.
ACM Trans. on Graphics (SIGGRAPH 2005)
24, 3 (2005), 1134–1141.G. Irving, J. Teran, and R. Fedkiw. 2004. Invertible Finite Elements for Robust Simulationof Large Deformation. In
Symp. on Computer Animation (SCA) . 131–140.A. Jacobson, I. Baran, J. Popović, and O. Sorkine. 2011. Bounded biharmonic weightsfor real-time deformation.
ACM Trans. on Graphics (TOG)
30, 4 (2011), 78.Petr Kadlecek, Alexandru-Eugen Ichim, Tiantian Liu, Jaroslav Krivanek, and LadislavKavan. 2016. Reconstructing Personalized Anatomical Models for Physics-basedBody Animation.
ACM Trans. Graph.
35, 6 (2016).Ladislav Kavan, Dan Gerszewski, Adam W. Bargteil, and Peter-Pike Sloan. 2011. Physics-Inspired Upsampling for Cloth Simulation in Games.
ACM Trans. Graph.
30, 4, Article93 (2011), 10 pages.Ali Emre Kavur, M. Alper Selver, Oˇguz Dicle, Mustafa Bariş, and N. Sinem Gezer. 2019.
CHAOS - Combined (CT-MR) Healthy Abdominal Organ Segmentation Challenge Data .https://doi.org/10.5281/zenodo.3362844Michael Kazhdan and Hugues Hoppe. 2013. Screened Poisson Surface Reconstruction.
ACM Trans. on Graphics (TOG)
32, 3, Article 29 (2013).Seunghwan Lee, Ri Yu, Jungnam Park, Mridul Aanjaneya, Eftychios Sifakis, and JeheeLee. 2018. Dexterous manipulation and control with volumetric muscles.
ACMTransactions on Graphics (SIGGRAPH 2018)
37, 4 (2018), 57:1–57:13.S. H. Lee, E. Sifakis, and D. Terzopoulos. 2009. Comprehensive Biomechanical Modelingand Simulation of the Upper Body.
ACM Trans. on Graphics
28, 4 (2009), 99:1–99:17.T. McInerney and D. Terzopoulos. 2008. Deformable Models. In
Handbook of MedicalImage Processing and Analysis (2nd Edition) , I. Bankman (Ed.). Chapter 8, 145–166.M. Müller and M. Gross. 2004. Interactive Virtual Materials. In
Proc. of Graphics Interface2004 . 239–246.M. Müller, B. Heidelberger, M. Teschner, and M. Gross. 2005. Meshless DeformationsBased on Shape Matching. In
Proc. of ACM SIGGRAPH 2005 . 471–478.G. Niculescu, J.L. Nosher, M.D. Schneider, and D.J. Foran. 2009. A deformable modelfor tracking tumors across consecutive imaging studies.
Int. J. of Computer-AssistedRadiology and Surgery
4, 4 (2009), 337–347.James F. O’Brien, Adam W. Bargteil, and Jessica K. Hodgins. 2002. Graphical Modelingand Animation of Ductile Fracture. In
Proceedings of ACM SIGGRAPH 2002 . 291–294.W. Press, S. Teukolsky, W. Vetterling, and B. Flannery. 2007.
Numerical recipes: The artof scientific computing (third ed.). Cambridge University Press, Cambridge, UK.T. Rhee, U. Neumann, J. Lewis, and K. S. Nayak. 2011. Scan-Based Volume AnimationDriven by Locally Adaptive Articulated Registrations.
IEEE Trans. on Visualizationand Computer Graphics
17, 3 (2011), 368–379.Shunsuke Saito, Zi-Ye Zhou, and Ladislav Kavan. 2015. Computational Bodybuild-ing: Anatomically-based Modeling of Human Bodies.
ACM Trans. on Graphics(SIGGRAPH 2015)
34, 4 (2015).J. Schmid, E. Gobbetti J. A. I. Guitián, and N. Magnenat-Thalmann. 2011. A GPUframework for parallel segmentation of volumetric images using discrete deformablemodels.
The Visual Computer
27 (2011), 85–95.Jérôme Schmid and Nadia Magnenat-Thalmann. 2008. MRI Bone Segmentation UsingDeformable Models and Shape Priors. In
Medical Image Computing and Computer-Assisted Intervention – MICCAI 2008 . 119–126.J. Schmid, A. Sandholm, F. Chung, D. Thalmann, H. Delingette, and N. Magnenat-Thalmann. 2009. Musculoskeletal Simulation Model Generation from MRI Data Setsand Motion Capture Data.
Recent Advances in the 3D Physiological Human (2009),3–19.Eftychios Sifakis, Igor Neverov, and Ronald Fedkiw. 2005. Automatic determination offacial muscle activations from sparse motion capture marker data.
ACM Trans. onGraphics (SIGGRAPH 2005)
24, 3 (Aug. 2005), 417–425.Breannan Smith, Fernando De Goes, and Theodore Kim. 2018. Stable Neo-HookeanFlesh Simulation.
ACM Trans. Graph.
37, 2 (2018), 12:1–12:15., Vol. 1, No. 1, Article . Publication date: August 2020. odeling of Personalized Anatomy using Plastic Strains • 17
Olga Sorkine and Marc Alexa. 2007. As-rigid-as-possible surface modeling. In
Symp. onGeometry Processing , Vol. 4. 109–116.O. Sorkine, D. Cohen-Or, Y. Lipman, M. Alexa, C. Rössl, and H-P Seidel. 2004. Laplaciansurface editing. In
Symp. on Geometry processing
ACM Trans. on Graphics (SIGGRAPH 2013)
32, 4 (2013),102:1–102:10.Robert W Sumner and Jovan Popović. 2004. Deformation transfer for triangle meshes.
ACM Trans. on Graphics (SIGGRAPH 2004)
23, 3 (2004), 399–405.J. Sylvester. 1884. Sur l’equations en matrices p x = x q.
C. R. Acad. Sci. Paris.
99, 2(1884), 67–71, 115–116.G. Székely, A. Kelemen, C. Brechbühler, and G. Gerig. 1996. Segmentation of 2-D and3-D objects from MRI volume data using constrained elastic deformations of flexibleFourier contour and surface models.
Medical Image Analysis
Proc. of ACM SIGGRAPH 1999 . 335–342.C. Twigg and Z. Kačić-Alesić. 2010. Point Cloud Glue: constraining simulations usingthe procrustes transform. In
Symp. on Computer Animation (SCA) . 45–54.Christopher D. Twigg and Zoran Kačić-Alesić. 2011. Optimization for Sag-Free Sim-ulations. In
Proc. of the 2011 ACM SIGGRAPH/Eurographics Symp. on ComputerAnimation
ACM Trans. on Graphics (SIGGRAPH2019)
38, 4 (2019).Bin Wang, Longhua Wu, KangKang Yin, Uri Ascher, Libin Liu, and Hui Huang. 2015b.Deformation capture and modeling of soft objects.
ACM Transactions on Graphics(TOG) (SIGGRAPH 2015)
34, 4 (2015), 94.Yu Wang, Alec Jacobson, Jernej Barbič, and Ladislav Kavan. 2015a. Linear subspacedesign for real-time shape deformation.
ACM Transactions on Graphics (TOG)(SIGGRAPH 2015)
34, 4 (2015), 57.Max A. Woodbury. 1950. Inverting modified matrices.
Memorandum Rept. 42, StatisticalResearch Group, Princeton University, Princeton, NJ (1950), 4pp.Hao Zhang. 2004. Discrete combinatorial Laplacian operators for digital geometryprocessing. In
Proceedings of SIAM Conference on Geometric Design and Computing
A PLASTIC STRAIN LAPLACIAN AND ITS NULLSPACE
Let L sc ∈ R m × m denote the discrete mesh Laplacian for scalar fieldson mesh tetrahedra [Zhang 2004], L sc i , j = i = j , − i (cid:44) j , and i , j are adjacent , . (28)Given a plastic strain state s = [ s , s , s , s , s , s ] ∈ R m , where s i ∈ R m , define the plastic strain Laplacian Ls = [ L sc s , √ L sc s , √ L sc s , L sc s , √ L sc s , L sc s ] , (29)where the √ s , s , s controltwo entries in the symmetric matrix F p ∈ R × . Lemma:
Assume that the tet mesh M has a single connected com-ponent. Then, the nullspace of L is 6-dimensional and consists ofvectors ψ i = [ s , . . . , s ] where s j ∈ R m is all ones when j = i , andall zeros otherwise. Proof:
First, observe that L sc is symmetric positive semi-definitewith a single nullspace vector, namely the vector of all 1s. Thisfollows from the identity x T L sc x = (cid:213) i and j adjacent ( x i − x j ) , (30) i.e., x T L sc x = x i are the same.We have 0 = s T Ls = (cid:205) i = ξ i s i T L sc s i , where ξ i = √ i = , ,
5; and 1 otherwise. Because L sc is symmetric positive semi-definite, each s i must either be or a non-zero nullspace vectorof L sc , i.e., a vector of all 1s. A linearly independent orthonormalnullspace basis emerges when we have a vector for all 1s for exactlyone i . There are 6 such choices, giving the vectors ψ i ; we normalizethem by dividing with √ m . ■ B FIRST AND SECOND DERIVATIVES OF ELASTICENERGY WITH RESPECT TO PLASTIC STRAIN
For convenience, we denote F p , i as i − th entry of the vector vec ( F p ) ∈R . The first-order derivatives are ∂ E ∂ x i = V ∂ ψ ∂ F e : ∂ F e ∂ x i = V P : ∂ F e ∂ x i , (31) ∂ E ∂ F p , i = V ∂ ψ ∂ F p , i + ∂ V ∂ F p i ψ = V P : ∂ F e ∂ F p , i + ∂ V ∂ F p i ψ , where (32) ∂ F e ∂ x i = ∂ F ∂ x i F − p , ∂ F e ∂ F p , i = F ∂ F − p ∂ F p , i , and (33) ∂ V ∂ F p , i = ∂ (cid:12)(cid:12) F p (cid:12)(cid:12) ∂ F p i V . (34)Here, P is the first Piola-Kirchhoff stress tensor and ∂ F / ∂ x is aconstant matrix commonly used in the equations for FEM simulation.For the second-order derivatives, we first compute ∂ E / ∂ x . Thisis the tangent stiffness matrix in the FEM simulation under a fixed F p . It is computed as ∂ E ∂ x i ∂ x j = V ∂ F e ∂ x j T : ∂ P ∂ F e : ∂ F e ∂ x i . (35)Here, ∂ P / ∂ F e is a standard term in FEM nonlinear elastic simulation;it only depends on the strain-stress law (the material model). Next,we compute ∂ E /( ∂ x ∂ F p ) ,∂ E ∂ x j ∂ F p , i = ∂ V ∂ F p , i (cid:18) P : ∂ F e ∂ x i (cid:19) + (36) V ∂ F e ∂ F p , j T : ∂ P ∂ F e : ∂ F e ∂ x i + V P ∂ F e ∂ x i ∂ F p , j , where (37) ∂ F e ∂ x i ∂ F p , j = ∂ F ∂ x i ∂ F − p ∂ F p , j . (38)Finally, we have ∂ E ∂ F p , i ∂ F p , j = ∂ V ∂ F p , i ∂ F p , j ψ + (39) V (cid:18) P : ∂ F e ∂ F p , i ∂ F p , j + ∂ F e ∂ F p , j T : ∂ P ∂ F e : ∂ F e ∂ F p , i (cid:19) + (40) ∂ V ∂ F p , j ∂ ψ ∂ F p , i + ∂ V ∂ F p , i ∂ ψ ∂ F p , j , where (41) ∂ V ∂ F p , i ∂ F p , j = ∂ (cid:12)(cid:12) F p (cid:12)(cid:12) ∂ F p , i ∂ F p , j V , (42) ∂ F e ∂ F p , i ∂ F p , j = F ∂ F − p ∂ F p , i ∂ F p , j . (43) , Vol. 1, No. 1, Article . Publication date: August 2020. The quantities ψ , P and ∂ P / ∂ F e are determined by the chosen elasticmaterial model. After computing the above derivatives, there isstill a missing link between F p and s . Because we want to directlyoptimize s , we also need the derivatives of E ( F p ( s ) , x ) with respectto s . From Equation 2 we can see that F p is linearly dependent on s . Therefore, so we can define a matrix Y such that vec ( F p ) = Ys .Then all the derivatives can be easily transferred to derivation by s by multiplying with Y . C PROOF OF SINGULAR LEMMA
Statement (i) follows from well-known linear algebra facts R( A ) = N( A T ) ⊥ and dim (N( A )) + dim (R( A )) = p , and the symmetry of A . As per (ii), A maps R( A ) into itself, and no vector from R( A ) mapsto zero, hence the restriction of A to R( A ) is invertible, establishinga unique solution to Ax = b with the property that x ⊥ ψ i for all i = , . . . , k . This unique solution is the minimizer ofmin x x T Ax − b T x (44)s.t. ψ Ti x = i = , . . . , k . (45)When expressed using Lagrange multipliers, this gives Equation 15.Suppose x = n + r is another solution and n ∈ N( A ) and r ∈R( A ) . Then b = Ax = Ar and hence r is the unique solution fromEquation 15. The vector n can be an arbitrary nullspace vector,proving the last statement of (ii). As per (iii), suppose we have0 = B ( n + r ) = Ar + (cid:205) ki = λ i α i ( ψ Ti n ) ψ i . Observe that the first summandis in R( A ) and the second in N( A ) . Hence, B ( n + r ) can only be zeroif both summands are zero. Ar = r = . The secondsummand can only be zero if n ⊥ ψ i for each i , which impliesthat n = . Hence, B is invertible. The last statement of (iii) can beverified by expanding (cid:0) A + (cid:205) ki = α i ψ i ψ Ti (cid:1) (cid:0) x + (cid:205) ki = λ i α i ψ i (cid:1) . ■ D PROOF OF NULLSPACE LEMMA
Fig. 24.
Illustration of the nullspaceproof.
The original forces F i sum tozero. We have G i = RF i ; note that R isthe same for all tets. Therefore, the ro-tated forces G i also sum to zero. Hencethere is no change in the internal elasticforce under a rotation, i.e., infinitesimalrotations are in the nullspace of K . We are trying to prove that K ( F p ( s ) , x ) is 6-dimensionalfor any x that solves f e ( s , x ) = . First, if x is a solution, thentranslating all vertices of theobject by the same constant3-dimensional vector is alsoa solution. This means thatvector ψ i : = ( e i , e i , . . . , e i ) isin the nullspace of K , where e i ∈ R is the i -th standard ba-sis vector, for i = , , . Now,suppose we rotate the objectwith an infinitesimal rotation X (cid:55)→ X + e i × X . Observe thatfor general plastic strains s , the elastic forces in each in-dividual tet are not zero evenin the equilibrium x ; but thecontributions of elastic forces on a tet mesh vertex from all ad-jacent tets sum to zero. As we rotate the object, the forces con-tributed by adjacent tets to a specific tet mesh vertex rotate by the same rotation in each tet. Therefore, as these forces sum tozero, they continue to sum to zero even under the rotation (see Fig-ure 24). This means that the vector of infinitesimal displacements ψ + i : = [ e i × x , e i × x , . . . , e i × x n ] induced by the infinitesimal rota-tion is in the nullspace of K , for each i = , , . Here, x i are the com-ponents of x = [ x , x , . . . , x n ] . The vectors ψ i , i = , , , , , , form the nullspace of K . ■ Finally, we inform the reader that the nullspace of K ( F p ( s ) , x ) isonly 3-dimensional if x is not an elastic equilibrium. In this case,only translations are in the nullspace. Infinitesimal rotations arenot in the nullspace because under an infinitesimal rotation, thenon-zero elastic forces f e rotate, i.e., they do not remain the same.The assumption of x being the equilibrium shape is therefore crucial(and is satisfied in our method). E SECOND DERIVATIVE OF POLAR DECOMPOSITION