Steepest Descent Paths On Simplicial Meshes of Arbitrary Dimensions

Download as pdf or txt
Download as pdf or txt
You are on page 1of 12

Steepest descent paths on simplicial meshes of arbitrary dimensions

Mattia Natalia,b , Marco Attenea,∗, Giulio Ottonelloc


a CNR-IMATI, Italy
b Universityof Bergen, Norway
c University of Genova, Italy

Abstract
This paper introduces an algorithm to compute steepest descent paths on multivariate piecewise-linear functions on
Euclidean domains of arbitrary dimensions and topology. The domain of the function is required to be a finite PL-
manifold modeled by a simplicial complex. Given a starting point in such a domain, the resulting steepest descent
path is represented by a sequence of segments terminating at a local minimum. Existing approaches for two and three
dimensions define few ad-hoc procedures to calculate these segments within simplexes of dimension 1, 2 and 3. Un-
fortunately, in a dimension-independent setting this case-by-case approach is no longer applicable, and a generalized
theory and a corresponding algorithm must be designed. In this paper, the calculation is based on the derivation of the
analytical form of the hyperplane containing the simplex, independently of its dimension. Our prototype implemen-
tation demonstrates that the algorithm is efficient even for significantly complex domains.

1. Introduction direction of the maximum speed is the direction of the


negative gradient of F. When the function is expressed
Being able to efficiently calculate the local minima through a closed analytical form, computing the gradi-
of a real-valued function is of clear importance for nu- ent becomes a matter of symbolic calculus and thus is
merous application contexts. In some domains, such as relatively easy. In contrast, when the function is de-
physics or chemistry, it is not rare to deal with prob- scribed through a finite set of samples, one may need to
lems where the function F to be minimized represents employ more sophisticated techniques.
the status of a system, and its local minima represent In this paper, we review the methods proposed so
particular configurations for which the system is stable far to compute steepest-descent paths on piecewise lin-
or in equilibrium. It is thus interesting to know how ear surface meshes, and generalize them within a uni-
such a system evolves towards a stable state when left fied approach that works on simplicial complexes of
unconstrained from an initial configuration. Such an any dimension. Furthermore, we provide a novel so-
evolution can be represented within the domain of F by lution to resolve the cases of indifferent equilibrium that
a sequence of configurations that form a path connect- may arise when tracking the steepest descent path. Fi-
ing the initial state with the final stable state. Normally nally, we show how our solution can be exploited in the
the system cannot jump from a configuration to another field of Material Sciences to study the crystallization of
without continuity (though there are exceptions, e.g. in molten substances.
quantum mechanics) and thus the path is continuous as Note that generalizing existing algorithms is not as
well. easy as it could appear at a first sight. On triangle
For several natural phenomena, a given non-stable meshes, for example, there are three possibilities: (1)
state tends to change towards equilibrium with maxi- the path passes through triangles, or (2) it may run
mum speed. Such a phenomenon can be modeled by along edges or (3) it may cross vertices. Thus, algo-
a real-valued function F that associates an energy with rithms treating these models can enumerate a few pos-
each state of the system and, if F is differentiable, the sibilities and handle each of them as a particular case.
Furthermore, the treatment of special cases of indif-
∗ CorrespondingAuthor
ferent equilibrium can also be handled by enumerating
Email address: marco.attene@ge.imati.cnr.it (Marco the possibilities. Conversely, in a generalized approach
Attene) such as the one introduced in this article, a dimension-
Preprint submitted to Computer & Graphics May 23, 2013
independent solution must be provided which can treat 1.1.2. Piecewise linear functions on simplicial meshes
all the cases based on a common approach. To the best When the domain is a two-dimensional simplicial
of our knowledge no solution to this problem has been mesh, several methods exist to calculate steepest de-
proposed so far, while there are application contexts scent paths. A big family of algorithms uses descent
(e.g. Material Sciences) that call for a solution. (and ascent) paths to bound the cells of the so-called
Morse-Smale complexes [3, 2, 1]. Within this family,
1.1. Related work the methods described in [4, 5, 6] force the steepest de-
scent paths to be sequences of mesh edges. Specifically,
The steepest descent method, also known as gradi-
these approaches start at a given vertex, and at each step
ent descent method, is an algorithm to compute local
simply choose the neighboring vertex associated with
minima of real-valued functions of n variables. While
the smallest function value. Though these algorithms
running, it constructs a piecewise-linear approximation
are extremely efficient, they produce only approximate
of the steepest descent path connecting an initial point
solutions. Conversely, in [7, 8] the paths are computed
to a local minimum.
based on [1] and are free to cross the triangles, therefore
produce much more precise results. These algorithms
1.1.1. Steepest descent direction for differentiable func- have applications in both Computer Graphics [9, 10]
tions and Visualization [11].
Formally, if f : Rn −→ R is a real function and Clearly, if the function is piecewise-linearly defined
∂f ∂f ∂f
∇ f (x) = ( ∂x , , . . . , ∂x
1 ∂x2 n
) is its gradient, then the steep- on a simplicial mesh, the gradient is generally unde-
est descent direction for x0 ∈ Rn is given by the vec- fined for points on edges and for the vertices, so it may
tor d = −∇ f (x0 ). When starting from an initial point be necessary to provide an estimate for these particular
x0 , the steepest descent method computes a new point cases. If pi is the position of the ith vertex of a closed
x1 := x0 − α∇ f (x0 ). Then, in a second step it com- and manifold triangle mesh M, then the gradient of the
putes another point x2 := x1 − α∇ f (x1 ), and so on. In function f : M −→ R can be estimated as [12]:
general, the algorithm computes a sequence of points
x0 . . . xk that is expected to converge to a local mini-
X
∇ f (pi ) := [ f (p j ) − f (pi )] w j , (1)
mum. The positive parameter α might be kept constant j∈N(i)
or can change from step to step. Small values of α pre-
where the w j s are proper weights that do not depend
vent the algorithm to oscillate around minima, but too
on f , and N(i) is the set of vertices connected to i by an
small values might significantly slow down the entire
edge.
process.
In contrast, the gradient ∇ f can be computed exactly
Near minima, the steepest descent method tends to
for points within triangles. Furthermore, ∇ f does not
go slowly and in some cases to have an erratic path. For
change across the points within a given triangle. So, on
these reasons, other minimization algorithms have been
a triangle t with vertices (pi , p j , pk ) and unit normal N,
proposed: notable examples are the conjugate gradient
the gradient is the solution ∇ f |t of the 3×3 linear system
that uses conjugate directions instead of the local gra-
[13]:
dient to take into account the previously-chosen direc-
tions, and the Newton-Raphson that exploits informa- 
 p j − pi 
  
 f (p j ) − f (pi ) 
tion of the second derivative to locate the minimum of  pk − p j  ∇ f |t =  f (pk ) − f (p j )  (2)
the function.
   
N 0
With small modifications, the steepest descent
method can be adapted to the case of piecewise-linear Therefore, alternatively to Eqn. 1, the steepest direc-
meshes [1] approximating differentiable functions. In tion at the vertex position pi can be estimated by consid-
this case, there is no need to choose the parameter α be- ering its incident triangle’s gradient having largest mag-
cause each segment of the path is bounded by the size nitude.
of mesh elements. As a consequence, the aforemen- The most detailed description of an algorithm that
tioned drawbacks are no longer an issue. Unfortunately, computes the steepest paths is given for the 2-
though there are some methods for calculation on two dimensional case in [1]. Therein, the path can go along
([1]) and three ([2, 3]) dimensional meshes, a complete edges or pass through faces of the triangulation. Three
adaptation is missing for the cases where differentiable cases occur for a starting point p of the path: (1) p is in
functions are defined on n-dimensional simplicial do- the interior of a triangle - the steepest direction is unique
mains. and orthogonal to the level lines (eg. Eqn. 2); (2) p is
2
in the interior of an edge - one or two locally steepest L is a subcomplex of K if L is a complex and L ⊂ K.
directions are possible; (3) p is a vertex - there may be For A ∈ K, the (closed) star of A in K, star(A, K), is the
as many locally steepest directions as the number of its subcomplex of K made of all simplices of K having A
incident triangles. as a face plus all their faces. If A ∈ K, then the link of
The case of 3-dimensional manifold domains is A in K, link(A, K), is the set of simplices in star(A, K)
treated in [2, 3], where smooth manifolds are approx- whose intersection with A is empty.
imated with piecewise linear simplicial complexes that A geometric realization |Ak | of a simplex Ak in the
linearly interpolate samples of the function f . Both in Euclidean space Rn , n ≥ k, can be obtained by defining
[3] and in [2] ascending/descending arcs are defined a bijection between the vertices of Ak and a set of k + 1
as piecewise linear curves (1-manifolds) between sad- affinely independent points p0 , p1 , . . . , pk of Rn , so that
|Ak | = {(t0 p0 + t1 p1 + . . . + tk pk ) ∈ Rn | ti ≥ 0, i ti = 1}.
P
dles and maxima/minima lying along edges of the input
mesh. Thus, |Ak | is the convex hull of p0 , . . . , pk and is said to
be an Euclidean simplex.
2. Definitions and problem statement The underlying space |K| of K is the union ∪A∈K |A|
of the geometric realization of its simplices. An under-
In this section basic definitions are introduced to sup- lying space |K| of a k-dimensional simplicial complex
port a formal description of the algorithm. The follow- K is unambiguously defined by a bijection between all
ing definitions are adapted from [14] and [15]. the vertices of K and a corresponding set of points of
Rn , n ≥ k, such that the image of the vertices of each
2.1. Simplicial complexes simplex consists of affinely independent points.
A k-dimensional simplex, or k-simplex, Ak is a set
V = {v0 , . . . , vk } of k + 1 objects called vertices, together 2.2. Problem statement
with the set of real-valued functions α : V → R satis- Herewith, some novel definitions are introduced to
fying vi ∈V α(vi ) = 1 and α(vi ) ≥ 0. A function α is
P
better explain the problem and the proposed algorithm.
called a point of Ak . The values α(v0 ), . . . , α(vk ) are the In particular, we need to formally define sets of com-
barycentric coordinates of the point α. binatorial entities that, after a geometric realization,
A (proper) face B of Ak , denoted B < Ak , is a simplex assume a special relationship with a given Euclidean
determined by the (proper) subset W ⊂ V, whose points point.
β : W → R are identified with the points α : V → R Let |K| be the underlying space of K, and let x be a
such that α(vi ) = β(vi ) if vi ∈ W and α(v j ) = 0 if v j ∈ point in it. We say that the influence set of x in K is
V − W. If B is a face of A, then A is said to be incident the set of all the simplices whose geometric realization
at B. contains x. That is:
A finite simplicial complex K is a finite set of sim-
plices such that: iset(x, K) = {σi ∈ K : x ∈ |σi |} (3)
i) if A ∈ K and B < A, then B ∈ K; Furthermore, we say that the boundary set of x in K
ii) if A, B ∈ K, then A ∩ B is either empty or it is a face is the set of all the faces of simplices in iset(x, K) whose
of both A and B. geometric realization does not contain x. That is:

From now on, we shall omit the term finite. bset(x, K) = {σi : σi < σ j , σ j ∈ iset(x, K), x < |σi |}
The dimension of K is the dimension of the largest (4)
dimensional simplex belonging to K. A simplicial com- For example, if we denote with t the geometric real-
plex of dimension n is homogeneous if it is made of ization of a 2-simplex σ ∈ K, t is a triangle, the influ-
n-simplices and their faces. Within an n-dimensional ence set of its barycenter is σ whereas its boundary set
simplicial complex, a k-simplex is said to be maximal if is made of the three vertices and the three edges being
k = n. faces of σ. Notice that the influence set of the geomet-
The boundary ∂A of a simplex A is the complex made ric realization of a vertex |v| coincides with its open star
of the proper faces of A. The boundary ∂K of a ho- (i.e. star(v) \ link(v)), and the boundary set of |v| is ex-
mogeneous n-dimensional simplicial complex K is the actly link(v). It is worth mentioning, however, that we
(n−1)-complex obtained as the sum mod 2 of the (n−1)- cannot simply replace iset and bset with the more stan-
dimensional simplices of the boundary ∂A of each of the dard star and link. Though a relation actually exists,
n-simplices A ∈ K plus their faces. in fact, iset and bset map arbitrary Euclidean points to
3
3. Algorithm description

In the easiest case, the algorithm starts with a given


point p0 within a (maximal) n-dimensional simplex σ0 ,
computes the direction of steepest descent −∇0 and
shoots a ray from p0 in the direction of −∇0 . The ray
(a) (b) intersects an n − 1-dimensional simplex which is a face
of both σ0 and an opposite maximal simplex σ1 . Let p1
be the point of intersection and −∇1 be the direction of
steepest descent calculated within σ1 . Once again, the
algorithm shoots a ray from p1 in the direction of −∇1 ,
intersects another n − 1-dimensional simplex, and so on,
as long as ∆i is positive.
Actually, this is just a particularly easy case. In
practice lower-dimensional simplices (i.e. of dimension
(c) (d)
< n − 1) can be intercepted along the path, can contain
Figure 1: Possible positions of a starting point in a two dimensional the starting point p0 , or can even contain whole parts of
complex. In (a), iset(p, K) = {σ1 } and just one direction must be the path. Thus, a procedure is necessary for the calcu-
evaluated to proceed to the next point q. When p is on an edge, the lation of the gradient vector within both maximal sim-
ith ray may (b), or may not (c), intersect a face of σi belonging to
plices (see section 4.1) and lower-dimensional simplices
bset(p, K) (note that in (c) each of the two arrows points outwards wrt
the corresponding σi ). If p is a vertex (d), the cardinality of iset(p, K) (see section 4.2). Furthermore, when the point pi lies
increases and several directions must be considered. on a simplex of dimension < n − 1, the gradient must
be computed for all its incident simplices of any dimen-
sion, and the path proceeds on the simplex for which ∆i
is maximized.
subcomplexes, whereas the star and link map simplexes Figure 1 illustrates the various possibilities for a start-
to subcomplexes, and a conversion makes sense only ing point p when K is two-dimensional. In (a), p lies
when the Euclidean point being mapped coincides with within a maximal simplex, so it is sufficient to calculate
the geometric realization of a vertex, as in the afore- the gradient within the simplex to determine the next
mentioned example. point q of the path, which might lie either on a ver-
tex or on an edge of the simplex. In (b), p is on a 1-
Let K be a finite n-dimensional simplicial complex, dimensional simplex, so it is necessary to compute the
VK be the set of its vertices, and r : VK → Rn be a map- gradient for the simplex itself and for both its incident
ping function defining an underlying space |K|r which triangles, and choose the one that maximizes ∆i ; in this
is an n-manifold with boundary (i.e. |K|r is a PL n- case ∆i is maximum for one of the incident triangles. In
manifold [14] with boundary). Also, let us consider a contrast, in (c) p is on a valley line, so ∆i is maximum
function f ∗ which associates a real value with each ver- on the 1-dimensional simplex. Note that in this case the
tex of K, and a corresponding map f : |K|r → R such ray shot from p towards the negative gradient of one of
that, if σ is a vertex x = r(σ) =⇒ f (x) = f ∗ (σ), the incident triangles would not intercept any edge of
whereas if σ is a higher-dimensional simplex and x ∈ the same triangle. Finally, in (d) p is on a 0-simplex, so
|σ| then the value of f (x) is defined by linearly interpo- all its incident simplices must be considered.
lating the value of f at its vertices.
3.1. Construction pipeline
Having said that, our problem can be summarized as
follows: given an initial point p0 ∈ |K|r , determine a Observe that the assumption that |K| is a manifold
sequence of points p0 . . . pk , pi ∈ |K|r , such that pi+1 ∈ with boundary guarantees that |bset(|v|, K)| is a topolog-
|bset(pi , K)|r maximizes the quantity ∆i = f (p i )− f (pi+1 ) ical sphere if v is an internal vertex, whereas it is a topo-
kpi+1 −pi k2
and ∆i > 0. logical semi-sphere if v is on the boundary. Therefore,
if v is internal, each ray emanating from it must neces-
Thus, the objective of the algorithm is not just the sarily intersect a point of |bset(|v|, K)|.
determination of the local minimum, but the determina- Let K be an n-dimensional mesh, with its geomet-
tion of the whole path that reaches that minimum while ric realization |K| that is an n-dimensional PL-manifold
following the direction of steepest descent. with boundary in Rn . Also, let p0 ∈ |K| be the starting
4
point of the steepest descent path we are looking for. Algorithm 1 Steepest descent path construction.
In the first step, the algorithm determines the set of sim- Require: A simplicial complex K, the value of f ∗ at
plices in the influence set of p0 , that is, iset(p0 , K). Note each vertex of K, and a point p0 on |K|.
that the cardinality of this set depends on the dimen- Ensure: A sequence p0 . . . pk , pi ∈ |K|, s.t. pi+1 ∈
sion of the minimal simplex whose realization contains |bset(pi , K)| maximizes ∆i and ∆i > 0.
p0 . For example, if p0 corresponds to the realization of
a vertex |v|, then all the simplices incident at v are in 1: i=0
iset(p0 , K). On the other extreme, if p0 belongs to the 2: ∆max = −∞
interior of the realization of a maximal simplex, then 3: compute iset(pi , K)
only that simplex will be in iset(p0 , K). In any case, 4: for each simplex σij in iset(pi , K) do
for each simplex σi in iset(p0 , K) we derive the gradi- 5: compute (∇ f )σ j
i
j
ent vector (∇ f )σi as described in section 4, and compute 6: pi+1 := intersection between |bset(pi , K)| and the
the intersection pi1 between |bset(p0 , K)| and the ray em- ray from pi in the direction of −(∇ f )σ j
i
anating from p0 in the direction of −(∇ f )σi . If such an 7:
j
if (pi+1 j
exists AND pi+1 is on a face of σij ) then
intersection does not exist, it means that p0 is on the j j
8: compute ∆i = ( f (pi ) − f (pi+1 ))/||pi+1 − pi ||2
boundary of K and the computed direction would move
9: if (∆i > ∆max ) then
the path outside the domain; thus, we declare −(∇ f )σi
10: ∆max = ∆i
to be an “invalid” direction. In general, a direction j
11: pi+1 = pi+1
−(∇ f )σi is declared to be “valid” if and only if the inter-
12: end if
section pi1 exists and belongs to a face of σi . For all the
13: end if
“valid” directions computed within iset(p0 , K), we com-
14: end for
pute the quantity ∆i = ( f (p0 ) − f (pi1 ))/||pi1 − p0 ||2 and
15: if ( ∆max > 0 ) then
select the point p1j that maximizes it. If such a maximum 16: i++
quantity is negative, we have reached a local minimum. 17: goto 2:
Otherwise, p1 = p1j becomes the new starting point for 18: end if
our path computation, and the algorithm iterates.
A pseudo-code describing the method is given in al-
gorithm 1. the path (in this order) and the algorithm iterates (Fig.
2(c)).
3.1.1. Flat simplices
This strategy to deal with flat simplices is based on
If the maximum ∆i is zero, it means that we have the assumption that f ∗ is typically computed by sam-
reached a “flat” simplex having the same function value pling another continuous function g using a regular pat-
associated with all its vertices. In the case of flat sim- tern. In cases such as those in Figures 3 and 12, for
plices, we have adopted a strategy based on the follow- example, it is easy to find flat simplexes approximating
ing principle: if a simplex σF is flat but some of its parts of g which have an unambiguous steepness. This
neighbors are not, then such a flatness may be a conse- is even more evident when the samples are taken on con-
quence of the discretization of the domain of f , so we tours: in Figure 6, for example, there are numerous flat
may consider it as a sort of artefact. In contrast, if both triangles that actually cut descending crests and valley
σF and all its neighbors are flat, then we assume that the lines. Note that in all these cases, isolated flat triangles
function is well represented and is actually flat in that can be there just because of the regularity of the pat-
region. Therefore, when the path reaches a flat simplex tern used, and thus our approach to cross them gives the
σF (Fig. 2(a)), we consider the union of the influence expected results.
sets of all its (realized) vertices, that we conveniently
denote as iset(|σF |, K) (Fig. 2(b)). For each simplex
σi ∈ iset(|σF |, K), σi , σF and σi ≮ σF , we check the 4. Gradient evaluation
flatness. If all the σi s are flat, then the algorithm termi-
nates. Otherwise, for each non-flat σi we compute the In order to implement the proposed algorithm, we
next point pF as the barycenter of the maximal common need a procedure to derive the gradient vector of f at
face of σi and σF , then compute the successive point each point in its domain |K|r . A point p ∈ |K|r can ei-
pF+1 and evaluate ∆i based on it. If the maximum ∆i ther be the image of a vertex or belong to the image
among all the σi s is negative the algorithm terminates, of a higher-dimensional simplex. We first consider the
otherwise the corresponding pF and pF+1 are added to case where p belongs to the image of a maximal sim-
5
(a) (b) (c)

Figure 2: Proceeding on a flat region: (a) a flat simplex σF and (b) its corresponding iset(|σF |, K); (c) pF is the edge’s midpoint (i.e. its barycenter)
and pF+1 maximizes ∆i .

plex (i.e. if K is n-dimensional, a simplex is maximal if piecewise-linear, the resulting n-dimensional Euclidean
its dimension is n). Afterwards, we generalize the result simplex spans a hyperplane Πσ of Rn+1 .
to lower-dimensional simplices through dimensionality To derive the analytic form of Πσ we exploit the
reduction. above determinant formulation by instantiating the
points Ai with the enriched vertices of |σ|, that is:
4.1. Evaluation in maximal simplices
f (x1 , . . . , xn )

x1 · · · xn 1
Our aim is to find the gradient vector of f at all 1
p1 · · · pn 1
f (p1 , . . . , pn )
1 1
1

the points within an Euclidean simplex |σ| defined by .
.. .. .. .. .. = 0.
n + 1 vertices, let them be p1 = (p11 , . . . , p1n ), . . . , pn+1 = . . . .
, . . . ,
n+1
1 , . . . , pn ). In this section we derive the analytic
(pn+1 n+1 p1 · · · pn+1
n f (p n+1
1 p n+1
n ) 1
form, inside |σ|, of the function f explicitly defined only
By expanding this determinant, for example with re-
at the vertices and assumed to be linear within the sim-
spect to the first row, and making it equal to zero, we
plex. Such an analytical form is a hyperplane, thus the
obtain:
gradient vector is constant and can be easily derived by
applying the definition: λ1 x1 + . . . + λn xn − λn+1 f (x1 , . . . , xn ) + λ0 = 0,
∂f ∂f ∂f λi ∈ R, which can be rewritten as
∇ f (x) = ( , ,..., ).
∂x1 ∂x2 ∂xn
λ1 λn λ0
f (x1 , . . . , xn ) = x1 + . . . + xn + .
To summarize, f : |σ| −→ R is explicitly defined just λn+1 λn+1 λn+1
at the vertices and linear interpolation is assumed for The above equation is the analytic form of f inside the
the other (internal) points. Then, the analytic form of f Euclidean simplex defined by points p1 , . . . , pn+1 , thus
within the simplex is computed (Section 4.1.1) and the its gradient vector is
n gradient vector coordinates are defined as the first n
coefficients of such an analytic form. ∂f ∂f λ1 λn
∇ f (x) = ( ,..., )=( ,..., ).
∂x1 ∂xn λn+1 λn+1
4.1.1. Hyperplane by n points
From an algorithmic point of view, the entire process
Let A1 , . . . , An be n affinely independent points of the requires the computation of n + 1 determinants of (n +
n-dimensional Euclidean space (eg. the vertices of an 1) × (n + 1) matrices. If we denote
Euclidean simplex). The implicit equation of the hyper-
plane spanned by A1 , . . . , An can be obtained by expand- f (p11 , . . . , p1n )
 1
 p1 · · · p1n 1 

M =  ... .. .. .. ..  ,



ing the determinant: 
 . . . . 
pn+1 n+1 n+1
, . . . , n+1 
· · · p f (p p ) 1
x1 · · · xn 1 1 n 1 n
← A1 → 1 thus λi is the determinant of the matrix obtained by

.. .. = 0,

. .

eliminating the i-th column from M.
← An → 1
Example: We assume n = 3 and we want to find
We observe that the value of f can be appended as the gradient vector of function f inside a simplex de-
an (n + 1)th coordinate to each vertex of |σ|. Since f is fined by the four points p1 = ( 12 , 32 , 1), p2 = (1, 1, 2),
6
p3 = (2, 2, 0) and p4 = (0, 0, 0), with respective values data structure with adjacencies [16] enriched with func-
f (p1 ) = 10, f (p2 ) = 2, f (p3 ) = 3, f (p4 ) = 0. Thus the tion values at vertices. Namely, if the mesh is n-
hyperplane equation (in R4 , where unknows are x, y, z dimensional and made of NV vertices and NS maximal
and w = f (x, y, z)) is defined by simplices, our data structure explicitly encodes:
• an array G of NV n-uples of coordinates represent-

x y z f (x, y, z) 1
1 3 1 10 1 ing the vertex positions along with an array F of
2 2
1 = 0. NV corresponding function values;

1 1 2 2

2 2 0 3 1
• an array S of NS (n + 1)-uples of indexes in G rep-
0 0 0 0 1

resenting maximal simplices along with an array O
that becomes of NS (n + 1)-uples of indexes in S representing the
corresponding adjacent maximal simplices;
−30x + 36y + z − 4 f (x, y, z) = 0.
• an array M of NV indexes in S representing, for
Therefore each vertex, one of the incident maximal simplices.
15 1 Using this structure, we can extract all the geometrical
f (x, y, z) = − x + 9y + z and topological information required for the execution
2 4
in optimal time. The algorithm has been implemented
and in C++.
15 1
∇ f (x) = (− , 9, ).
2 4
4.2. Evaluation in lower-dimensional simplices
When the simplex is not maximal, the approach de-
scribed in Section 4.1 cannot be applied “as it is” be-
cause the matrix would be rectangular. Thus, our ap-
proach is based on a projection of the simplex onto a
lower-dimensional space so that it becomes maximal. (a) (b)
In such a new space we can derive the gradient as de-
Figure 3: Path computed on a known function for verification. The
scribed in Section 4.1, and the resulting vector can be
surface is the image of f (x, y) = sin(x) + cos(y). The path starts from
transformed back to the original high-dimensional space (0, 0) in (a) and from (1, 1) in (b).
by inverting the projection matrix.
Thus, let |σ| be defined by m vertices p1 =
(p1 , . . . , p1n ), . . ., pm = (pm
1
1 , . . . , pn ), with m < n + 1.
m

We observe that the set of m − 1 vectors g1 = p2 − p1 ,


g2 = p3 − p1 , . . ., gm−1 = pm − p1 constitutes a basis
of the (m − 1)-dimensional Euclidean space where σ is
maximal. So, we translate |σ| so that its first vertex p1
coincides with the origin of Rn . Then, we consider the
(m−1)×(n) matrix G made of the gi s and orthonormalize
it (eg. through Gram-Schmidt). By right-multiplying
the orthonormal matrix G⊥ by each vertex of |σ| we ob-
tain an (m − 1)-dimensional Euclidean simplex in Rm−1 .
After having derived the (m − 1)-dimensional gradi-
ent vector of f as described in Section 4.1, we right-
multiply GT⊥ by ∇ f to transform it back to the original
space Rn . Figure 4: Steepest descent path computed on a large model of
the Kilimanjaro volcano (1 million triangles, Model courtesy of
Aim@Shape’s repository). On volcanos, these paths can be used to
predict possible lava flows.
5. Implementation and results

Since the algorithm is dimension-independent, in our Though we have tested our prototype on meshes of
implementation we have employed an extended indexed various dimensions, for n ≥ 4 visualizing the results is
7
tricky and not very informative. Thus, with the excep-
tion of figure 12, we herewith show some of our results
for domains of dimensions two and three only. Note
that in some two-dimensional cases the value of f may
represent a height (eg. Figures 4, 6 and 3), whereas in
some other cases the function may have a completely
different meaning (eg. a temperature, such as in Figure
5). Nevertheless, to better conceive the results through
our illustrations, for any two-dimensional mesh we have (a) (b)
used the value of f as a third coordinate, thus embed-
ding the 2-dimensional complex in R3 , even if f does Figure 7: A cubical domain with an associated function f (x, y, z) = z.
In (a) the cube is aligned with the coordinate system, and the path tra-
not represent a height. In contrast, for most illustrations
verses only lower-dimensional simplices. In (b) the cube is rotated: as
regarding three-dimensional domains (Figures 7, 8, 9 expected, the path descends vertically through the solid until it reaches
and 10) we have used the height (i.e. the z coordinate) the boundary, then it follows the steepest direction on the boundary of
to define the value of f , with the exception of Figure 11 the domain.
where the function values are provided aside and repre-
sent a potential.

Figure 5: The so-called ”liquidus surface” of a ternary chemical sys-


tem. The local minimum reached by the path represents a stable state Figure 8: Two points of view of the same path computed on a tetrahe-
of the material. dral mesh. Model courtesy of Aim@Shape’s repository.

6320 @ 1.87GHz and 4096MByte RAM. In table 1 and


table 2 execution times are reported for meshes of di-
mension two and three respectively. To provide the
reader an idea of the complexity of the input, the ta-
bles indicate both the number NV of vertices and the
number NS of maximal simplices. However, the num-
ber of operations (and hence the execution time) per-
formed by the algorithm depends mostly on the number
of simplices traversed by the path.

Figure 6: Steepest descent path on 2-dimensional simplicial complex


representing a terrain.

5.1. Computational times


The size of the experimental data sets shown in this
paper ranges from a few hundred to more than a million
simplices. Experiments have been performed on a stan- Figure 9: Paths computed starting from two different points of the
dard desktop PC equipped with a Intel(R) Core(TM)2 same tetrahedral mesh. Model courtesy of Aim@Shape’s repository.

8
Model NV NS time
CAS (fig.5) 3777 7416 0.027
Quark (fig.3) 3969 7688 0.019
Vobbia (fig.6) 7745 15312 0.051
Kilimanjaro (fig.4) 490001 977204 0.78

Table 1: Execution time for the models of dimension two made of NV


vertices and NS triangles. Time is expressed in seconds and does not
include input/output operations.

Model NV NS time
Figure 10: A path computed starting from an internal point of a dense
Cube (fig.7(a)) 182 436 0.042 tetrahedral mesh. As expected, the path goes vertically as long as it
Rotated Cube (fig.7(b)) 103 239 0.026 can, while it correctly follows the boundary when necessary. Model
Aorta (fig.8) 9529 35551 0.219 courtesy of Pierre Alliez.
Horse (fig.9) 837 2152 0.076
Dino (fig.10) 11979 61527 0.231
Bucky (fig.11) 262144 1250235 0.456 h is the sampling step used, we may consider the follow-
ing cases: if h is too large some simplices might actu-
Table 2: Execution time for the models of dimension three made of ally cancel small local minima, and the resulting path
NV vertices and NS tetrahedra. Time is expressed in seconds and does might be longer than expected; if h is sufficiently small
not include input/output operations.
but x0 is closer than h to the border of a maximal cell
of the stable Morse complex of f (e.g. Fig 3(a)), the
5.2. Robustness, stability and accuracy path may take a wrong direction and fall into a different
minimum; in all the other cases the path ends within a
Instability is intrinsic at each point of the domain
maximum distance h from the ideal minimum.
where the gradient vanishes, and small perturbations
may lead to dramatically different gradient directions. Furthermore, when the samples of f are corrupted by
In particular, if the point is a local maximum of a differ- noise, the function itself should be smoothed to pre-
entiable function, any direction in the domain is equally vent the path to be broken at a spurious local mini-
eligible to track a steepest-descent path. The most con- mum. While advanced techniques exist for two dimen-
servative approach would be to compute all the possible sions [18, 19], the case of higher dimensions has been
paths, but unfortunately they might be an infinite num- mostly left uninvestigated. Finally, the computation of
ber. If f is a Morse function [17] algorithms exist to the hyperplane of a nearly degenerate simplex is ill-
compute the so-called descending Morse complex; from conditioned, and also in this case existing work to re-
such a structure one can easily understand which are the move such undesired degeneracies is mostly dedicated
points that can be possibly part of a descent path starting to two and three dimensions [20].
at a local maximum. Unfortunately, these algorithms When the function is actually piecewise linear (i.e.
are known only for domains of dimension two and three. it is not an approximation) our algorithm provides an
For this reason, for now we just ignore the problem and exact result. This happens, for example, in the case
let the algorithm provide one of the possible solutions. of a linear program whose constraints define a convex
From the computational point of view, the discretiza- polytope. In this case one may also use Dantzig’s sim-
tion might be a further issue. In several cases, a smooth plex algorithm while being guaranteed to reach the same
function f is only partially defined through a finite sam- correct minima. Differently from the simplex method,
pling of its domain which is linearly interpolated within however, our algorithm is free to cut through higher-
a simplicial mesh. To be conservative, also in this case dimensional simplexes and can work on nonconvex do-
we may need to produce infinite solutions where the mains. On convex domains, passing through edges only
path crosses a flat simplex. However, as mentioned in (as done by the simplex algorithm) is sufficient to reach
section 3, depending on the surroundings we may inter- the final exact minimum, though the path is not neces-
pret such a flatness as an artefact and let the algorithm sarily exact. In contrast, the same constraint on noncon-
provide one of the possible solutions. vex domains might cause the path to reach completely
It is also worth speculating about the accuracy of our wrong results.
algorithm when it is used to process a discretized ver- In case of 2- and 3-dimensional meshes, existing
sion of a Morse function f . If x0 is the starting point and methods can be applied to compute steepest paths.
9
These paths are sometimes called integral curves and a vector specifying the direction of movement when ly-
they are mainly used to generate Morse-Smale com- ing in it, and the magnitude of the vector corresponds
plexes. Amongst these existing algorithms, some [7, 8] to the steepness in that direction. In this setting, for
allow paths to cut through maximal simplices, as we each simplex σ the value of ∆i is just the magnitude
do, and perform in similar time. In these cases our al- of the vector vσ associated with σ. If the input vectors
gorithm achieves exactly the same results while being form a nongradient vector field [21] some of the paths
more general. Conversely, other methods (e.g. [4], [5], may be closed loops, thus the visited simplices should
[6]) are more efficient than ours, but only return an ap- be marked to avoid the algorithm to infinitely loop.
proximation of the steepest path since they build it as a
collection of edges (1-simplices). For this reason, these
methods often require an initial mesh refinement to re-
duce inaccuracy. 5.4. Applications in material research

Material research (ceramurgy, glass technology, slags


industry), Geology (petrology, geochemistry) and en-
gineering try to understand the characteristics of com-
plex materials by studying the interrelationship of com-
position, microstructure and process conditions repre-
sented in phase diagrams, whose computation is one of
the main objectives of computational thermodynamics.
Figure 11: A complex tetrahedral mesh (1.25 million tets) discretizing In these fields of research the absolute amount of pure
a potential f around a carbon molecule with a path reaching a local components constituting a complex material is not very
minimum of f . Model courtesy of Aim@Shape’s repository. important, while it is much more interesting to know
their relative amount. For example, pseudowollastonite
is always composed of 50% CaO and 50% S iO2 in mo-
5.3. Generalization lar proportions independently of the bulk amount of the
The algorithm described in this paper has numerous substance. A convenient way to describe the domain of
practical applications, including the prediction of lava all the possible relative compositions is by using a sim-
flow paths (Fig. 4) and the computation of stable states plex, where each vertex corresponds to one of the pure
reachable by cooling molten substances (Fig. 5). In components. Thus, such a simplex is called the compo-
some cases, however, the knowledge of the function sitional simplex. For any molten substance, the temper-
alone is not sufficient to accurately predict the evolu- ature of incipient crystallization changes depending on
tion of a system, and further variables may be necessary the relative quantity of pure components. The function
to provide a more comprehensive and precise model of that associates the incipient crystallization temperature
the phenomenon. For example, when trying to predict a with each point of the compositional simplex is called
lava flow, one should take into account that the flowing the liquidus. Since the dimension of the simplex is arbi-
material has an inertia. Thus, to better predict the path, trary (i.e. it depends on the number of pure components)
we should correct the steepest descent direction by us- the liquidus is a hypersurface and, as such, can be repre-
ing an additional inertial term. As a further example, sented piecewise-linearly through a simplicial complex
in material sciences one can predict the result of a crys- [22].
tallization process by following descent lines on the so- By starting at a given point in the compositional sim-
called liquidus hypersurface (see Section 5.4), but such plex (i.e. a relative composition) it is possible to track
lines are not necessary locally steepest at all their points, so-called descent lines representing the crystallization
because their direction depends on additional variables. induced by loss of heat toward the exterior or by com-
We may adapt our algorithm by providing a general- pression. At each point, the direction of such lines is
ized approach to determine the direction to be followed dictated by the lowering of freezing point equation [23].
at all the points. Specifically, the gradient field can be Thus, by using this equation we may assign the correct
computed in advance (once for each simplex) and possi- direction to each simplex of the liquidus, and then ex-
bly corrected by considering additional terms. Thus, the ploit the algorithm described here to predict the compo-
input to the adapted version of the algorithm is a sim- sition that will be generated when the process reaches
plicial complex where each simplex is associated with its stability point.
10
require a substantially different algorithm.

Acknowledgements. This work is partly supported by


the Italian MIUR-PRIN Project N. 2009B3SAFK
(Topology of phase diagrams and lines of descent) and
t t by the Petromaks program of the Norwegian Research
Council through the Geoillustrator project (N. 200512).
Thanks are due to the SMG members at IMATI for help-
z
z
y
x
ful discussions.

References
[1] H. Edelsbrunner, J. Harer, A. Zomorodian, Hierarchical morse
z complexes for piecewise linear 2-manifolds, in: SCG ’01: Pro-
t ceedings of the seventeenth annual symposium on Computa-
tional geometry, ACM, New York, NY, USA, 2001, pp. 70–79.
[2] H. Edelsbrunner, J. Harer, V. Natarajan, V. Pascucci, Morse-
y y
x smale complexes for piecewise linear 3-manifolds, in: SCG ’03:
x
Proceedings of the nineteenth annual symposium on Computa-
tional geometry, ACM, New York, NY, USA, 2003, pp. 361–
Figure 12: An example showing a path computed within a four- 370.
dimensional domain. The function is defined as f (x, y, z, t) = (y − [3] V. Natarajan, V. Pascucci, Volumetric data analysis using morse-
x2 )2 + (x − 1) + z + t. The domain is a unitary hypercube [−0.5, 0.5] × smale complexes, in: SMI ’05: Proceedings of the International
[−0.5, 0.5] × [−0.5, 0.5] × [−0.5, 0.5], and the images depict four axis- Conference on Shape Modeling and Applications 2005, IEEE
aligned 3D projections of the result. To better convey the result, edges Computer Society, Washington, DC, USA, 2005, pp. 322–327.
of the simplicial complex are not shown. [4] S. Takahashi, T. Ikeda, Y. Shinagawa, T. L. Kunii, M. Ueda,
Algorithms for extracting correct critical points and construct-
ing topological graphs from discrete geographic elevation data,
6. Conclusions and future work Computer Graphics Forum 14 (3) (1995) 181–192.
[5] C. Bajaj, D. Schikore, Topology preserving data simplification
with error bounds, Comput. Graph. 22 (1) (1998) 3–12.
This research shows that it is possible to efficiently [6] P. Bremer, V. Pascucci, A practical approach to two-dimensional
compute steepest descent paths on piecewise-linear scalar topology, in: Topology-Based Methods in Visualization,
functions defined on Euclidean domains of any dimen- Springer, Berlin Heidelberg, Germany, 2007, pp. 151–169.
[7] P.-T. Bremer, H. Edelsbrunner, B. Hamann, V. Pascucci, A
sion. In practice, such functions are defined through a
multi-resolution data structure for two-dimensional morse func-
finite sampling and the resulting discretization may lead tions, in: VIS 03: Proceedings of the IEEE Visualization 2003,
to flat simplices for which we have proposed an appro- IEEE Computer Society Press, Los Alamitos, CA, 2003, pp.
priate treatment. In some cases, a vector field may be 139–146.
[8] V. Pascucci, Topology diagrams of scalar fields in scientific vi-
available instead of the function, and often such a field sualization, in: Topological Data Structures for Surfaces, John
determines the path’s direction at all the points; in these Wiley and Sons London, U.K., 2004, pp. 121–129.
cases, we have shown that our algorithm can be eas- [9] S. Dong, P.-T. Bremer, M. Garland, V. Pascucci, J. C. Hart,
ily adapted to compute the paths according to the input Spectral surface quadrangulation, ACM Trans. Graph. 25 (3)
(2006) 1057–1066.
field, and have presented a practical application in the [10] X. Ni, M. Garland, J. C. Hart, Fair morse functions for extracting
area of Material Sciences. the topological structure of a surface mesh (2004) 613–622.
As it is, however, the algorithm presented here has [11] J. L. Helman, L. Hesselink, Representation and display of vector
some limitations that constitute the subject of our future field topology in fluid flow data sets, Computer 22 (8) (1989)
27–36.
research. First, when the steepness is the same in more [12] S. Biasotti, G. Patanè, M. Spagnuolo, B. Falcidieno, Analysis
than one direction at a given point our algorithm sim- and comparison of real functions on triangulated surfaces, Curve
ply chooses one of the possibilities, whereas it would and Surface Fitting Modern Methods in Mathematics (2007) 41–
be certainly useful to have a more comprehensive result 50.
[13] S. Biasotti, G. Patanè, M. Spagnuolo, B. Falcidieno, G. Bare-
including possible branches and joins: the possibility quet, Shape approximation by differential properties of scalar
to have even infinite alternatives (eg. when the starting functions, Computers and Graphics 34 (3) (2010) 252–262.
point is within a flat simplex) makes this problem rather [14] M. Attene, D. Giorgi, M. Ferri, B. Falcidieno, On converting
sets of tetrahedra to combinatorial and PL manifolds, Computer-
challenging. An even more challenging direction for fu- Aided Geometric Design 26 (8) (2009) 850–864.
ture research is to consider higher-degree (i.e. non lin- [15] L. C. Glaser, Geometrical Combinatorial Topology, Vol. 1, Van
ear) functions to interpolate the samples, but this might Nostrand Reinhold, 450 West 33rd Street, New York, 1970.

11
[16] L. De Floriani, A. Hui, Data structures for simplicial complexes:
an analysis and a comparison, in: SGP ’05: Proceedings of the
third Eurographics symposium on Geometry processing, Euro-
graphics Association, Aire-la-Ville, Switzerland, Switzerland,
2005.
[17] S. Biasotti, L. De Floriani, B. Falcidieno, P. Frosini, D. Giorgi,
C. Landi, L. Papaleo, M. Spagnuolo, Describing shapes by
geometrical-topological properties of real functions, ACM
Comp. Surveys 40 (4) (2008) 12:1–12:87.
[18] H. Edelsbrunner, D. Morozov, V. Pascucci, Persistence-sensitive
simplification of functions on 2-manifolds, in: Procs. of Sympo-
sium on Computational Geometry, 2006, pp. 127–134.
[19] G. Patane, B. Falcidieno, Computing smooth approximations
of scalar functions with constraints, Computers and Graphics
33 (3) (2009) 399–413.
[20] M. Attene, A lightweight approach to repairing digitized poly-
gon meshes, The Visual Computer 26 (11) (2010) 1393–1406.
[21] G. Chen, K. Mischaikow, R. S. Laramee, P. Pilarczyk, E. Zhang,
Vector field editing and periodic orbit extraction using morse
decomposition, IEEE Trans. on Visualization and Computer
Graphics 13 (4) (2007) 769–786.
[22] M. Attene, G. Ottonello, Computational geometry meets mate-
rial science, ERCIM News 84 (2011) 43–44.
[23] G. Ottonello, Principles of Geochemistry, Columbia University
Press, N.Y., 1997.

12

You might also like