Steepest Descent Paths On Simplicial Meshes of Arbitrary Dimensions
Steepest Descent Paths On Simplicial Meshes of Arbitrary Dimensions
Steepest Descent Paths On Simplicial Meshes of Arbitrary Dimensions
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.
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
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
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.
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
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
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