Using simulated annealing to obtain good nodal approximations of deformable bodies
Using Simulated Annealing to Obtain Good Nodal Approximations of Deformable Bodies Oliver Deussen and Leif Kobbelt and Peter T¨u cke
Institute for Operating and Dialog Systems,University of Karlsruhe
76128Karlsruhe,Germany
Net:[oliver|kobbelt|tuecke]@informatik.uni-karlsruhe.de
WWW:https://www.docsj.com/doc/d05272997.html,rmatik.uni-karlsruhe.de/[?oliver|?kobbelt|?tuecke]
Abstract.In this paper we present a method to obtain good approximations of de-
formable bodies with spring/mass systems.An iterative algorithm based on voronoi
diagrams is used to get a good mass distribution.The elastic properties of the system
are optimized by simulated annealing.Results are shown,and some applications
are discussed.
Keywords:
Simulation,Spring/mass lattice,Modeling,Deformable Bodies,ComputerGraphics Real time simulation and interactive control of deformable bodies are of great interest to computer graphics and other?elds like pattern matching or medical visualization.
Spring/mass(or nodal)systems are often used to obtain ef?cient simulations of deformable objects if lower requirements for the simulation precision are given.
Terzopoulos and Waters[15,13]used them for modeling and simulating facial tissue.Simulation of woven clothes was done by Breen[4].In an early work,Miller [11]used nodal models to simulate utricular objects.
In medical research similar systems are used to simulate human organs for training tasks in endoscopic surgery[10,6],and for general surgery simulation and planning[5].
To reproduce a speci?c mechanic behaviour of a spring/mass system,parameters like spring constants can either be preset,or optimized iteratively.We show that within some general limits,such optimizations yield systems which approximate objects with given elastic properties quite well.
1Nodal Approximations
In general,a spring/mass system(SMS)is a collection of points(or particles)with a speci?c mass,that are connected by springs,dampers or plastic elements.The discrete Lagrangian equation of motion for a point at position and mass has the simple
form(cf.[15]):
¨˙1
where is a general damping coef?cient,and denotes an external force acting on the point.The forces are exerted by elastic,viscotic and plastic connectors between and neighbour points(:set of neighbours).
1.1Finite Element Formulation of Nodal Systems
The idea of approximating a deformable body with a system of springs is surprisingly old.In1868,Kirsch[9]derived the fundamental equations of elasticity from a system of springs approximating a small cube.His work was one of the bases of the?nite element method(FEM).
Today pin jointed trusses(assemblages of axially loaded elastic bar elements)are used in structural mechanics as idealizations for many problems.The FEM approach is widely used for such systems.
Those trusses can be seen as spring/mass systems.The only difference is that within trusses bars have a mass and joints are massless whereby in a SMS massless springs are connected by mass points that behave like spherical joints.In both cases,the mass distribution can be expressed by a mass matrix.
If only elastic and viscotic connectors are used,equation(1)can be written in matrix calculus(cf.[2])by:
¨˙2 In this equation,111is the vector of all point positions1,is the damping matrix and denotes the external forces acting on the nodes.is the systems stiffness matrix.This matrix describes the forces on nodes as a result of deformations and can be seen as a generalized spring constant.
In comparison to other systems of?nite elements,nodal systems do have the advan-tage of being diagonal,which is simple and convenient for mathematical handling (cf.[8]).In addition,equation(2)can be solved by local application of equation(1) to each node.This leads to more ef?cient algorithms(cf.[15]),especially if parallel computers are used.
Table1shows some computation rates(calculations per second).The local simula-tion method based on equation(1)is independent of the bandwidth of.FEM methods based on equation(2)do have low computation rates if is recalculated in every step. This has to be done to make the FEM method comparable with the local algorithm.If
is updated every eighth step,the computation rates for systems with a small bandwidth are close to that of the local method.
For small deformations(the geometric linear case)can be assumed constant.In this case the computation rates are high,but simulation precision decreases.
The calculation of equilibrium states is computational expensive because a linear equation system based on has to be solved.
Table1:Computation rates(calculations per second,SUN Sparc20) Bandwidth of3721
3737
810
1931
2543
2.57
3.1The Iterated Voronoi Approach
Our observation was,that systems with arbitrary given points do not lead to good results in simulating elastic behaviour of homogeneous bodies.To solve this problem we developed a method of moving the points slightlyin order to obtain a more homogeneous point positioning.
For each given point
,the corresponding V oronoi region is de?ned by 3
The set of
resulting regions for all is a tesselation of the plane,it can be computed
in
log time for given points in two-dimensional space (cf.[12]).Since the boundary of each is polygonal,we can easily intersect it with
.(a)(b)(c)(d)
Figure 1:Point positioning -a)random points;b)V oronoi regions;c)cutted regions with moved points;d)regions after 15iterations,together with point motion.
Next,the are moved to the center of gravity of the corresponding V oronoi region (Figure 1).This process can be iterated.During iteration,the standard deviation of the point to point distances reduces as well as the standard deviation of the V oronoi-region sizes (Figure 2).
01020
30Iterations 0.000.05
0.10s t a n d a r d d e v i a t i o n i n [m ]0102030
Iterations 0.000.010.02s t a n d a r d d e v i a t i o n i n [m 2]Figure 2:Standard deviation during iteration -a)point to point distances;b)sizes of V oronoi-regions.
Both results indicate a more homogeneous point positioning after some iterations.After the last iteration the weights are computed according to section 4(see below).
4
3.2Regular Grids
Besides random positioning of points,we used regular grids together with speci?c connection topologies.These con?gurations lead to good approximation results with the investigated loads(see below),but may have an anisotropic behaviour because of their regularity(like crystals do).
Figure3:a)regular grids used;b)Suitable topologies to connect the grids.
To combine regular grids with arbitrary surfaces,tesselation is necessary.In the context of our quantitative study in the following paragraphs,we omit such tesselations because their results are hard to interpret.
4Second Step:Calculation of the Mass Distribution
If the point positions are known,we have to approximate the mass distribution of by assigning masses:to the points.The approximation of by a set of tupels is considered to be good,if the mass moments up to order two coincide.These moments determine the linear and angular accelerations to external forces and torques.
The relevant mass moments of are computed by
24
Substituting by a?nite set of points corresponds to approximating this integral by the sum
25
1
Equation(5)can be considered as the application of a cubature formula for the integrand function over the domain,where the points are the sample points and are the weight coef?cients for the formula.Hence,the approximation of meets all relevant mass moments if and only if this cubature formula has at least quadratic polynomial precision,i.e.,if and only if
6 where11is the vector of mass coef?cients,000100
is the10-dimensional vector of exact mo-010001110101011200020002
5
ments of,and is the10-matrix with
11111212121
12222222222
.. ...
.
..
.
..
.
..
.
..
.
..
.
..
.
1222
The solution of the under determined system(6)with least Euclidian norm2is of the form with10[3].Hence,the mass coef?cients can be found by solving
7 Since is a V andermonde matrix,is regular if there is a collection of10points which does not lie on a polynomial surface0of degree less than or equal to two.
Capturing third or higher order moments of yields no signi?cant improvement of the approximation und thus the degrees of freedom for the choice of the are used in a way that2is minimal.In this case the masses are distributed as uniform as possible.
5Third Step:Optimizing Elasticity
Given points with positions and masses.The connections can be found as follows: If the points are distributed randomly,we use a Delaunay triangulation to obtain pairs of points to be connected by springs.Two points are joined if they have a common edge in the V oronoi diagram.This triangulation method is widely used in?nite element
con?gurations used,is a square plate.In situations(1)and(2),stretching loads are applied,shearing is done in situations(3)and(4).
6
5.1Getting Reference Values
For some simple bodies and special con?gurations it is possible to get an analytic solution of the resulting displacements due to external forces (cf.[2]).2
As a reference,we calculate the deformation of the square plate shown in Figure 4(1).In this situation a uniform distributed force is applied to one side,the opposite one is constrained in -direction,but not in -direction.(a)
Figure 5:a)Plate with uniform distributed load to two opposite sides;b)resulting deformation.This situation is equivalent to the application of a uniform distributed load to two opposite sides,whereby the body is constrained such that rigid body motions are permitted (Figure 5).
We assume a linear elastic deformation and a material that can be expressed by Poisson’s ratio and Y oung’s modulus .In section 2we gave the general formulation of the strain-stress relation (see equation (3)).This relation,now written in matrix calculus,has to be satis?ed:
1
2For more complicated bodies and situations more precise ?nite element methods are used to get reference values.
7
which satis?es equation(9).The stresses can be computed by application of Hooke’s law(equation(8)):
Integration according to equation(10)leads to the desired displacements
We set and to zero,the proof of correctness for doing this is omitted here.
5.2Initial SMS Con?guration
The initial spring constants of the SMS are set according to Young’s modulus().For a spring of resting length0,the spring constant is set heuristicially to04. The upper part of Figure8(a)shows the calculated displacements of the initial SMS, approximating the deformations of the given surface.They are drawn using dotted lines.
5.3Fast Calculation of Displacements
During the optimization process,the displacements of the SMS which are due to the applied forces have to be calculated very often.Within the equilibrium state of the system they are obtained by solving
11 which is the solution of equation(2)with vanishing¨and˙.The additional vector describes the bearing forces in?xed points.The equation system is solved on the base of band matrices,the necessary node enumeration is done with an optimization method in a preprocessing step.Because is changed in every optimization step,it can not be inverted in advance.
5.4Quality Criteria
For each point position of the SMS and each of the test con?gurations,the reference displacement is calculated according to paragraph5.1.The aforementioned FEM approach is used to compute the actual displacements during optimization.
One possibility to describe the quality of an SMS approximation is the standard deviation between actual and reference displacements of all its points:
:
5.5Optimization by Simulated Annealing
Simulated annealing(cf.[1,14])is used as optimization method.This algorithm was originally designed for discrete problems but was later adapted to continuum problems. The method seems to be adequate because a pure gradient-descent delivered bad results, and analytic methods were very inef?cient because of the high degree of freedom within the system.Another advantage is the simple algorithm(Figure6).
temperature:=initial
of
steps:=0;
k:=0;
repeat
Choose randomly an edge;
Change randomly spring constant of;
Calculate according to equation(12);
if()then begin the system is now a better approximation :=;
good steps+1;
end else begin
Choose random number01;
if(
steps>glimit)
end
Figure6:Simulated annealing algorithm.
A drawback of the method is the high sensitivity to the number of temperatures, the initial temperature,the reduction function and the number of good steps(glimit) allowed in each temperature.Although,we found generally applicable values for those parameters(e=number of edges):
Parameter temperature nlimit
V alue097099
The simulated annealing process converges if nearly no bad steps are done any more (due to temperature)and the improvement at each temperature is below some threshold. We used an over-exponential temperature reduction to get convergence.
Another possibility is to stop the algorithm after a speci?c number of temperatures(e.g. no temperatures=2*e)as done in Figure6.In this case,the time complexity of the algorithm is2.
9
6Extension to3-D
Figure12shows the in?uence of V oronoi iterations on homogeneity and optimization quality of a nodal approximation.
8Conclusion and Future Work
Spring mass systems are often used for real time simulation of complex deformable objects.Our method allows the generation of systems with de?nite mechanic behaviour. We found a way to obtain an adequate mass distribution and simulated homogeneous materials as well as inhomogeneities and anisotropies.Systems up to some hundred points were optimized successfully.For large systems the optimization is computational expensive but this has to be done only once.
Besides the full implementation of the3-D case,we are working on frequency analysis.This is a promising approach because natural frequencies give some more information about the dynamic behaviour of deformable systems.Natural frequencies can be obtained by solving the eigenvalue problem
013 The frequencies are derived by2.Once again we have the advantage that
is diagonal,what makes ef?cient calculation methods applicable.
Frequency analysis is to be inserted into our optimization process in order to obtain nodal systems with good equilibrium and natural frequency approximations.
In addition,natural frequencies can be used to?nd very ef?cient integration methods for speci?c systems and situations.No stepsize control is needed if the highest natural frequencies are integrated with respect to the sampling theorem.This can be used in virtual reality applications like surgery simulation or virtual sculpting,for which our optimization process as well is of high interest.
References
1.E.H.L.Aarts and P.J.M.van Laarhoven.Statistical cooling:A general approach to
combinatorial optimization problems.Philips Journal of Research,(40):193–226, 1985.
2.K.-J.Bathe and E.L.Wilson.Numerical methods in?nite element analysis.
Prentice-Hall,1976.
3.W.Boehm and H.Prautzsch.Numerical methods.AK Peters Wellesley,1993.
4.D.E.Breen,D.H.House,and P.H.Getto.A physicially-based particle model of
woven cloth.The Visual Computer,(8):264–277,1992.
5.H.Delingette,G.Subsol,S.Cotin,and J.Pignon.A craniofacial surgery simulation
testbed.Research report2199,Institute National de Recherche en Informatique et Automatique,1994.
6.O.Deussen and Chr.Kuhn.Echtzeitsimulation deformierbarer Objekte¨u ber
nodale Modelle.in:Th.Strothotte and P.Lorenz,Hrsg.,Proc.Integration von Bild,Modell und Text.ASIM Mitteilungen No.46,University of Magdeburg, 1995.
11
7.J.Eisley.Mechanics of elastic structures.Prentice Hall,1989.
8.Xiaochun Gao,Zhiying King,and Qixian Zhang.A hybrid beam element for
mathematical modelling of high-speed?exible linkages.Mech.Mach.Theory, 24(1):29–36,1989.
9.G.E.Kirsch.Die Fundamentalgleichungen der Theorie der Elastizit¨a t fester
K¨o rper,hergeleitet aus der Betrachtung eines Systems von Punkten,welche durch elastische Streben verbunden sind.VDI-Zeitschrift,1868.
10.U.G.K¨u hnapfel,B.Neisius,H.G.Krumm,and M.H¨u bner.CAD-Based simulation
and modelling for endoscopic surgery.Endoscopic Surgery,1993.
https://www.docsj.com/doc/d05272997.html,ler.The motion dynamics of snakes and https://www.docsj.com/doc/d05272997.html,puter Graphics,
22(4):169–173,1988.
12.F.P.Preparata and https://www.docsj.com/doc/d05272997.html,putational geometry.Springer-V erlag New
Y ork,175Fifth Avenue,New Y ork,New Y ork10010,USA,3Au?age,1990. 13.D.Terzopoulos and K.Waters.Analysis and synthesis of facial image sequences
using physical and anatomical models.IEEE Transactions on Pattern Analysis and Machine Intelligence,1993.
14.P.J.M.van Laarhoven and E.H.L.Aarts.Simulated annealing,theory and appli-
cations.Reidel Publishing,Dordrecht,1987.
15.K.Waters.A physical model of facial tissue and muscle articulation derived
from computer tomography data.Visualization in Biomedical Computing,pp.
574–583.SPIE,1992.
Figure8:Initial(?rst line)and optimized displacements of a system with25nodes and connec-tion structure1.Stiffnesses for horizontal,vertical and diagonal springs are set to the same value (1:6:18:28055)
12
Figure9:If all spring constants are optimized separate(which is to be done for inhomogeneous bodies)the system approximates the given deformations(1:3:15:2306)but has some irregularities within the grid.
Figure10:Modeling anisotropic behaviour with rectangular cells(1:3:15:105048).The system has different elasticities in and-direction,thus the deformations due to the test loads are different.
Figure11:Simulation of an inhomogeneous material.The Young’s modulus of the left third of the plate is twice as high as the right part.In the upper part of the?gure the initial approximation of the spring/mass system was omitted.The lower part shows the optimized system.
13
a)Initial given points
b)System after one iteration
c)System after?ve iterations
d)System after15iterations
e)System after50iterations
V oronoi iterations
before optimization
1:3:15
0.450.470.440.420.37
Figure12:Optimized displacements with different number of preceding V oronoi-iterations.The The starting quality raises more than the quality after the optimization with100temperatures.
14