C Energy
Resume:
Discrete Comput Geom (****) **: *** ***
A Monotonicity Property for Weighted Delaunay
Triangulations
David Glickenstein
Received: 19 August 2005 / Revised: 8 May 2007 /
Published online: 11 September 2007
Springer Science+Business Media, LLC 2007
Abstract Given a triangulation of points in the plane and a function on the points,
one may consider the Dirichlet energy, which is related to the Dirichlet energy of a
smooth function. In fact, the Dirichlet energy can be derived from a nite element
approximation. S. Rippa showed that the Dirichlet energy (which he refers to as the
roughness ) is minimized by the Delaunay triangulation by showing that each edge
ip which makes an edge Delaunay decreases the energy. In this paper, we introduce
a Dirichlet energy on a weighted triangulation which is a generalization of the energy
on unweighted triangulations and an analogue of the smooth Dirichlet energy on a
domain. We show that this Dirichlet energy has the property that each edge ip which
makes an edge weighted Delaunay decreases the energy. The proof is done by a direct
calculation, and so gives an alternate proof of Rippa s result.
Keywords Triangulations Weighted Delaunay triangulation Dirichlet energy
Laplacian Rippa
1 Introduction
In [22], Rippa proved a remarkable theorem about the optimality of Delaunay triangu-
lations. He considered a quantity he called the roughness R(f, T ), which is a function
of some values f on a collection of data points in the plane and a (two-dimensional)
triangulation T of the data points. The roughness is de ned as the Sobolev semi-norm
2 2
i i
R(f, T ) = + d xdy,
x y
Ti
i
D. Glickenstein University of Arizona, Tucson, AZ, USA
e-mail: ************@****.*******.***
652 Discrete Comput Geom (2007) 38: 651 664
where i is the linear interpolation of f over the triangle Ti in T and the sum is over
all triangles in the triangulation. One may consider changing the triangulation by
exchanging two triangles joined by an edge, forming a quadrilateral, by the triangles
obtained by switching the diagonal of the quadrilateral; this is called an edge ip
or a 2 2 bistellar ip. He showed that the roughness of a triangulation decreases
when an edge is ipped to make the edge Delaunay. Since every triangulation can be
transformed into a Delaunay triangulation by a sequence of edge ips, this implies
that the roughness is minimized by a Delaunay triangulation.
The roughness is equal to the following functional, which we call the Dirichlet
energy:
1
E(f, T ) = (cot ij + cot ij )(fj fi )2,
4
{i,j } T1
where f is a function on the vertices of the triangulation T, T1 are the edges in the
triangulation, and ij, ij are the two angles opposite the edge {i, j } (one of the terms
is zero if {i, j } is on the boundary). This formula is often called the cotan formula
and goes at least as far back as Duf n s paper on lumped networks in 1959 [7]. It
has been popularized in the last ten years or so by the work of Pinkall and Polthier
[20]. Using the notation of the Dirichlet energy, Rippa s theorem may be stated as
follows. The precise de nitions of Delaunay triangulation and bistellar ip are given
in Sect. 2.
Theorem 1 (S. Rippa [22]) Let T be a triangulation of a nite number of points
in R2 . Let T0 be the vertices of the triangulation and let f : T0 R be a function.
Suppose T is another triangulation of the same points which is gotten from T by
a 2 2 bistellar ip of an edge e (in particular, T0 = T0,) such that the hinge is
Delaunay after the ip. Then
E(f, T ) E(f, T ),
where E(f, T ) and E(f, T ) are the Dirichlet energies corresponding to T and T .
As a consequence, the minimum is attained when all edges are Delaunay (and hence
the triangulation is a Delaunay triangulation).
The key step in Rippa s proof is the exact calculation of E(f, T ) E(f, T ). See
also Powar s method of calculating this quantity [21].
In this paper we generalize Rippa s monotonicity lemma to the class of weighted
Delaunay triangulations, sometimes referred to as regular triangulations or coherent
triangulations in the literature (see, for instance, [1, 8, 10, Chapt. 7]). In Sect. 2 we
recall how a set of weights can give a notion of the geometric dual (or center) of a
triangle. Given an edge, the centers of the two triangles incident on the edge may be
connected to form a dual edge which has a signed length. This allows us to give the
following de nition of Dirichlet energy.
Discrete Comput Geom (2007) 38: 651 664 653
De nition 2 The Dirichlet energy, depending on a triangulation T, a function f on
the vertices of T, and a set of weights w on the vertices, is de ned as
{i, j }
1
E(f, T, w) = (fj fi )2, (1)
{i, j }
2
{i,j } T1
where {i, j } denotes the length of edge {i, j } and {i, j } denotes the (signed)
length of the edge dual to {i, j }.
Given this de nition, the main theorem is the following generalization of Rippa s
monotonicity lemma.
Theorem 3 Let (T, w) be a weighted triangulation of some points in R2 . Let T0 be
the vertices of the triangulation and let f : T0 R be a function. Suppose (T, w) is
another weighted triangulation which is gotten from (T, w) by a 2 2 bistellar ip
on edge e such that the hinge is weighted Delaunay after the ip. Then
E(f, T, w) E(f, T, w),
where E(f, T, w) and E(f, T, w) are the Dirichlet energies of f corresponding to
(T, w) and (T, w).
The proof we provide has the advantage of being a straightforward calculation, but
the disadvantage that this calculation is quite long and complicated. It is an interesting
observation that the Dirichlet energy de ned for weighted triangulations does not
appear to come from a piecewise linear interpolation except in the case of equal
weights, when its formula is the cotan formula. The formula is derived in parallel to
the classical Dirichlet energy instead of as an approximation to it.
The organization of the paper is as follows. We begin by introducing weighted
Delaunay triangulations and the notation used in the rest of the paper. We then state
and prove the ip monotonicity result. Finally we comment about what would be
needed for a global theorem similar to Rippa s.
2 Weighted Delaunay Triangulations
In this section we x notation and review the de nitions of weighted Delaunay trian-
gulations. We denote simplices with braces { } such as {i } for a vertex, {i, j } for an
edge, and {i, j, k } for a triangle. We shall often omit the braces for a vertex, denoting
it simply as i . The entire triangulation is denoted T = {T0, T1, T2 }, where T0 are the
vertices (0-dimensional simplices), T1 are the edges (1-dimensional simplices), and
T2 are the triangles (2-dimensional simplices). In practice we may specify T by only
specifying T2, where the other pieces can be seen as the vertices and edges in the
triangles listed.
Two triangles adjacent to a common edge is called a hinge, and always has the
form {{i, j, k }, {i, j, }}. A triangulation can be altered by bistellar ips, the most
important for our purposes being 2 2 bistellar ips, which we will call edge ips or
654 Discrete Comput Geom (2007) 38: 651 664
Fig. 1 Hinges T and T
differing by a bistellar ip. The
dual edges are also shown
just ips. A ip replaces the hinge {{i, j, k }, {i, j, }} by the hinge {{i, k, }, {j, k, }}.
Such a ip can be seen in Fig. 1.
We shall study the triangulations of a given set of points T0 in R2 . The lengths of
edges in the triangulation (gotten as the distance between the points in R2 ) is given
by a function
: T1 (0, ),
where we write ij as the length of edge {i, j } T1 . The triangulation has the con-
cepts of volume of a one-simplex (length) {i, j } = ij and two-simplex (area of
a Euclidean triangle) {i, j, k } . If a hinge is convex, then a ip makes sense and
changes the function, where the new length can be calculated as the length of the
other diagonal.
We shall now add the weights to the structure.
De nition 4 A weighted triangulation (T, w) is a triangulation T together with
weights
w : T0 R.
We now recall the de nition of a weighted Delaunay triangulation (see, for in-
stance, [1] or [8]). We rst de ne centers which give geometric structure to the
Poincar dual of the triangulation. We will only need the (signed) length of edges
{i, j } dual to edges {i, j }, although one can de ne signed volumes of all of the
Poincar duals of simplices (in any dimension) as is done in [13].
Let d(x, p) be the Euclidean distance between points p and x . De ne the power
distance
p : R2 R
by
p (x) = d(x, p)2 wp (2)
if p is a point weighted with wp . The center C({i, j }) of edge {i, j } is de ned to be
the point x on the line containing {i, j } where i (x) = j (x). The center C({i, j, k })
Discrete Comput Geom (2007) 38: 651 664 655
Fig. 2 A weighted triangle with
circles corresponding to the
weights at the vertices. Central
orthogonal circle and dual edges
are also shown
of triangle {i, j, k } is the point x where i (x) = j (x) = k (x). The center of a
triangle can be given a weight according to the formula
wC({i,j,k }) = d({i }, C({i, j }))2 + d(C({i, j }), C({i, j, k }))2 wi .
Note that this number is independent of permutation of the indices. In the sequel,
it will become important to de ne signed distances between centers, but the sign of
the distances is irrelevant in the formula for the weight of the center because of the
squares in the formula.
If the weights are positive, weighted triangulations can be interpreted geometri-
cally (see Fig. 2). The power is then a function which is zero on the circle centered at
p with radius wp, positive outside the circle, and negative inside the circle. Notice
that if p is a vertex of a simplex and c = C then c (p) = wp and p (c) = wc .
If the weight wc is positive, it is the square of the radius of the circle orthogonal to the
vertex circles. If the weight is negative, there is still an interpretation as a degenerate
orthogonal circle using the methods in [19, Sect. 40], which describes the space of
circles as a vector space with an inde nite inner product generalizing the notion of
angle between circles.
It will become important to de ne a signed distance between centers as follows.
The distance from a triangle center C({i, j, k }) to an edge center C({i, j }) may be de-
ned to be the positive distance from C({i, j, k }) to C({i, j }) if the center C({i, j, k })
is in the same half plane determined by {i, j } as {i, j, k } is and to be negative one
times that distance if the center is in the opposite half plane than the triangle {i, j, k }.
We denote by hij,k this signed distance from C({i, j, k }) to C({i, j }), referring to the
fact that hij,k is like a height of {i, j } inside {i, j, k }. One can calculate hij,k explic-
itly as follows. Use a Euclidean motion to move {i, j, k } so that i is at the origin and
{i, j } is along the positive x -axis. Then hij,k is the y -component of C({i, j, k }). The
components of C({i, j, k }) can be found using the linear algebra of circles described
656 Discrete Comput Geom (2007) 38: 651 664
in [19, Sect. 40], which results in the following:
wj
wi wk
ij
hij,k = cot kij + cot j ik + cot ij k (3),
2 2
2 2Aij k
ij ij
where ij = {i, j } is the length, Aij k = {i, j, k } is the area, and ij k is the interior
angle at vertex i in triangle {i, j, k }. You can see the altitudes corresponding to hij,k
in Figs. 1 and 2. The dual to edge {1, 2} in Fig. 1 has negative length. Note that the
signed distance is de ned so that ij hij,k + i k hik,j + j k hj k,i = 2Aij k .
Recall De nition 2 for the Dirichlet energy. The (signed) length of the edge dual
to {i, j }, denoted as {i, j } , is equal to hij,k + hij, if {i, j } is an interior edge and
hij,k if {i, j } is on the boundary. Notice that the dual edge {i, j } is orthogonal to
the corresponding edge {i, j }. For a function f : T0 R on a weighted triangulation
of a hinge T = {{i, j, k }, {i, j, }}, the Dirichlet energy is equal to
hij,k + hij, hik,j hj k,i
E(f, T, w) = (fj fi )2 + (fk fi )2 + (fj fk )2
2 ij 2 ik 2 jk
hi hj,i,j
+ (f fi )2 + (fj f )2 .
2 2j
i
The lengths {i, j } being positive is equivalent to the following classical de n-
ition of weighted Delaunay.
De nition 5 An edge {i, j } adjacent to the two triangles {i, j, k } and {i, j, } is
weighted Delaunay if C({i,j,k }) > w and C({i,j, }) (k) > wk, where C({i, j, k })
is the center of {i, j, k } and similarly for {i, j, }. If the weights at the vertices are all
equal to zero, a weighted Delaunay edge is said to be Delaunay.
Sometimes we will instead say that the hinge is weighted Delaunay. A hinge is
Delaunay if and only if it satis es the local empty circumcircle property: the circle
circumscribing {i, j, k } does not contain . This is simply the interpretation of the
de nition when the weights are equal to zero. It is easy to see that an edge {i, j } con-
tained in a hinge {{i, j, k }, {i, j, }} is weighted Delaunay if and only if {i, j } 0.
The proof is essentially the same as the corresponding theorem for Delaunay trian-
gulations, which goes back at least to Rivin [23]. Details may be found in [12] or
simply proved by looking at when the length of the dual edge is equal to zero. The
formula for hij,k makes it easy to see that the condition for being weighted Delaunay
is unchanged by a weight scaling of the type {wi wi + c}i T0, where c is some
constant independent of i, since hij,k is unaffected by such a deformation.
The formula (1) for the Dirichlet energy is motivated as follows. We wish to write
the Dirichlet energy as
1
E(f, T, w) = fi fi Ai,
2
Discrete Comput Geom (2007) 38: 651 664 657
where is a Laplacian operator and Ai = i is the area of the dual to the vertex i .
We can write a general form for a Laplacian in terms of integration by parts
f
f dA = dS,
n
U U
where f is the normal derivative and dS is the measure on the boundary. In our case
n
we take U = i and this formula takes the form
{i, j }
fi Ai = (fj fi )
{i, j }
j
if we choose the dual structure properly. This leads to the formula (1). This is the
derivation of the Laplacian using Discrete Exterior Calculus as introduced in [14].
When the centers of triangles lie inside the triangles, it is possible to show that so-
lutions to the discrete Dirichlet problem associated to these Laplacians will actually
approximate solutions of the Dirichlet problem [6]. Special cases of Laplacians with
this form have also been studied, for instance, in [2, 3, 11, 12, 17, 18, 20]. These
Laplacians are all Laplacians on the graph determined by the one-skeleton of the tri-
angulation although the coef cients may be negative. In [7], Duf n interprets this
type of Laplacian in terms of an electrical network on the one-skeleton of the trian-
gulation. Each wire (edge) {i, j } has a resistor with resistance {i, j } {i, j } . This
make sense since the resistance of a wire should be proportional to its length and
inversely proportional to its cross-sectional area. See [4] for general theory of graph
Laplacians.
3 Monotonicity for Weighted Delaunay Triangulations
We now prove the main theorem, Theorem 3, for the Dirichlet energy on weighted
triangulations.
The proof depends on the following generalization of Rippa s key lemma [22,
Lemma 2.2] (see also [21]). Refer to Fig. 1 for a picture of the hinges.
Lemma 6 Let T ={{1, 2, 3}, {1, 2, 4}} and T = {{1, 3, 4}, {2, 3, 4}} be two hinges
differing by a ip along {1, 2}. Then
E(f, T, w) E(f, T, w) = (fT (c) fT (c))2 A2
1234,
where
2(r3 r4 r1 r2 )A1234 + w1 A234 + w2 A134 w3 A124 w4 A123
= (4),
8A123 A134 A234 A124
Aij k is the area of {i, j, k },
A1234 = A123 + A124 = A134 + A234 (5)
658 Discrete Comput Geom (2007) 38: 651 664
is the area of the hinge, c is the intersection of the diagonals, ri is the distance be-
tween c and vertex i, and fT and fT are the piecewise linear interpolations of f
with respect to the different triangulations. One can write
r1 r2
fT (c) = f2 + f1,
12 12
r3 r4
fT (c) = f4 + f3 .
34 34
The proof is somewhat involved, although straightforward. We use a proof which
is more direct than the ones given by Rippa [22] and Powar [21] for the case of
Delaunay triangulations.
Proof Recall the de nition of hij,k from (3), which we think of as the height of the
triangle {i, j, C({i, j, k For any function f, we can compute
4
1
E(f, T, w) E(f, T, w) = aij fi fj,
2
i,j =1
where
h12,3 h12,4 h13,2 h13,4
a12 = + a13 =,,
12 12 13 13
h14,2 h14,3 h23,1 h23,4
a14 = a23 =,,
14 14 23 23
h24,1 h24,3 h34,1 h34,2
a24 = a34 =,,
24 24 34 34
and aii = j =i aij (where we have symmetrized aij = aj i ). We now wish to factor
the coef cients.
We can easily gure out ri in terms of areas in the following way. Let vi be
the point in the plane representing {i }. We see that c = v1 + r112 (v2 v1 ) = v3 +
r3
(v4 v3 ). By taking the cross product with v2 v1 or v4 v3 we nd that
13
12 A134 34 A123
r1 = and r3 =,
A1234 A1234
(recall the de nition of A1234 from (5)) . Similarly,
12 A234 34 A124
r2 = and r4 = .
A1234 A1234
Thus
r3 r4 r1 r2
fT (c) fT (c) = f4 + f3 f2 f1
34 34 12 12
1
= (A123 f4 + A124 f3 A134 f2 A234 f1 ).
A1234
Discrete Comput Geom (2007) 38: 651 664 659
Also useful will be the calculation
1
r3 r4 r1 r2 = 2 (
2 2
(6)
34 A123 A124 12 A234 A134 ).
A1234
There are essentially two different types of coef cients to consider. We need only
consider a12 and a13 since the others are similar. Let ij k be the angle at vertex i in
triangle {i, j, k }. Consider a12 .
h12,3 h12,4
a12 = +
12 12
1 w1 w2 w3
= cot 312 + 2 cot 213 + 2 cot 123
2 4A123
2 12 2 12
1 w1 w2 w4
+ cot 412 + 2 cot 214 + 2 cot 124
2 4A124
2 12 2 12
1 w1
= (cot 312 + cot 412 ) + 2 (cot 213 + cot 214 )
2 2 12
w2 w3 w4
+ 2 (cot 123 + cot 124 ) (7)
.
4A123 4A124
2 12
Let be the angle at c in the triangle {1, 3, c}. We shall use the fact that in any triangle
{i, j, k } we have ij = i k cos ij k + j k cos j ik to compute the parts.
13 23 cos 312 14 24 cos 412
cot 312 + cot 412 = +
2A123 2A124
2 2
12 13 cos 123 12 14 cos 124
13
= + 14
2A123 2A124
2 2
sin 314
13
= + cot
14
2A123 2A124 sin sin 123
sin 413
+ + cot
sin sin 124
2 2
1 sin 314 sin 413
13
= + +
14
2A123 2A124 sin sin 123 sin 124
2 2
1 12 A134 A1234
13
= +
14
2A123 2A124 sin 34 A123 A124
+ 2 A123
2A 2A
13 124 12 134
14
=
2A123 A124
+2 + 2 A2 2 A2
2 2 A2 2A
12 134 A234
13-14-13-124-** 123 12 134
= +
2A1234 2A123 A124 A1234 2A123 A124 A1234
660 Discrete Comput Geom (2007) 38: 651 664
since
sin 314 = cos 123 sin + sin 123 cos
and
sin 413 = cos 124 sin sin 124 cos .
Furthermore,
+
2 2 2 2 2 2
13 A124 14 A123 12 A134
1
= + sin2 123 sin2 ( 123 + 124 ))
222 2
12 13 14 (sin 124
4
12 2 2
= (sin 123 sin 124 cos 134 )
2 12 13 14
= 2A123 A124 13 14 cos 134
since
sin2 A + sin2 B sin2 (A + B) = 2 sin A sin B cos(A + B).
Thus, using (6), we have
+ 2 13 14 cos 134 )
2 2 2A
( 12 134 A234
13 14
cot 312 + cot 412 =
2A1234 2A123 A124 A1234
2 A134 A234
2A
34 123 A124 12
=
2A1234 A123 A124
A1234
= (r3 r4 r1 r2 ).
A123 A124
For the other parts,
cos 213 cos 214
cot 213 + cot 214 = +
sin 213 sin 214
2A
sin 234 234
= = 12
sin 213 sin 214 2A123 A124
and
2A
12 134
cot 123 + cot 124 = .
2A123 A124
Thus (7) implies that
A234 A134
a12 = (2A1234 (r3 r4 r1 r2 ) + w1 A234
4A123 A134 A234 A124
+ w2 A134 w3 A124 w4 A123 ) = 2A234 A134
using the de nition of in (4).
Discrete Comput Geom (2007) 38: 651 664 661
Now consider a13 . We can compute
h13,2 h13,4
a13 =
13 13
1 w1 w3 w2
= cot 213 + 2 cot 312 + 2 cot 123
2 4A123
2 13 2 13
1 w1 w3 w4
cot 413 + 2 cot 314 + 2 cot 134
2 4A134
2 13 2 13
1 w1
= (cot 213 cot 413 ) + 2 (cot 312 cot 314 )
2 2 13
w3 w2 w4
+ 2 (cot 123 cot 134 ) + .
4A123 4A134
2 13
We see that
sin 324 sin 124
cot 213 cot 413 =
sin 213 sin sin 413 sin
2 A123 A124
2A
12 134 A234 34
=
2A1234 A123 A134
sin 324 = cos sin 213 + sin cos 213 sin 124 =
since and similarly
cos sin 413 + sin cos 413 . We also get
cos 312 sin 314 cos 314 sin 312
cot 312 cot 314 =
sin 312 sin 314
2A
sin 324 234
= 13
=
sin 312 sin 314 2A123 A134
and
2A
13 124
cot 123 cot 134 = .
2A123 A134
And so, using (6) and (4),
A234 A124
a13 = (2A1234 (r3 r4 r1 r2 ) + w1 A234
4A123 A134 A234 A124
w3 A124 + w2 A134 w4 A123 ) = 2A234 A124 .
A similar argument gives the other coef cients. Then we see, for instance, that
a11 = a12 a13 a14
= 2( A234 A134 + A234 A124 + A234 A123 ) = 2A2
234
with similar expressions for a22, a33, and a44 . Finally, we get that
E(f, T, w) E(f, T, w) = (A123 f4 + A124 f3 A234 f1 A134 f2 )2,
which is equivalent to the lemma.
662 Discrete Comput Geom (2007) 38: 651 664
Now we can prove the theorem.
Proof of Theorem 3 Since the coef cient a12 = {1122} and a34 = {3344} , we see
{, {,, } , }
that a12
Contact this candidate