MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

RE: DelaunayTriangulation

  • To: mathgroup at
  • Subject: [mg53210] RE: [mg53081] DelaunayTriangulation
  • From: "Melko, Mike" <Mike.Melko at>
  • Date: Thu, 30 Dec 2004 01:38:15 -0500 (EST)
  • Sender: owner-wri-mathgroup at

-----Original Message-----
From: peteraptaker [mailto:psa at] 
To: mathgroup at
Subject: [mg53210] [mg53081] DelaunayTriangulation

I hope do "DelaunayTriangulation" and operate on associated triangle
elements. I was hoping that DiscreteMath`ComputationalGeometry` would do
the job.

Question 1: is there a simple way of identifying the actual triangler
(rather than just the connections)

Question 2: If my desired outer bound is not convex how can I get my
derired triagulation and boundary nodes (no longer a convex hull I

Please excuse by finite element and Dealauney ignorance:



Mike Melko wrote:


I'm writing a package that uses "DelauneyTriangulation" for certain
and have spent some time thinking about this. Here is my attempt to
your questions. Let me say that I have found this a bit confusing, and
to articulate. Hopefully, the following explanation will be of some help
(to me
at least). Let me add the caveat that I'm a novice at computational

Let me start with Question 2 because the answer's short:

Question 2: If my desired outer bound is not convex how can I get my
triangulation and boundary nodes (no longer a convex hull I assume).

Answer: I actually forgot about this point, and I thank you for calling
to my attention. I believe the answer is that you must be able to
your boundary points somehow, and delete the unwanted adjacencies
procedurally from the output of "DelauneyTriangulation".

You might be able to do this by augmenting your data with "outside"
and then deleting those points and their adjacencies from the output of


Now for question 1:

Question 1: Is there a simple way of identifying the actual triangular
(rather than just the connections)

Answer: First, let me point out that I think the justification for the
structure used in the output of "DelauneyTriangulation" is that it
allows for
efficient searching and alteration of vertex adjacencies.

I'm not aware that DiscreteMath`ComputationalGeometry` provides any
to do what you want to do. I wrote a function that does something like
you're asking, because I also want to get at the actual triangles. Below
is a
simplified (and untested -- I don't have access to Mathematica at the
version of a function called "Triangulate". What follows after the
function is
a "dissection of the code. Unfortunately, it's rather long-winded.

Let me underscore that I'm not sure the following is correct. I intend
to test
it in the next couple of days to make sure. Hopefully it is a good


(*	Func:	Triangulate

	Desc:	Produces triangulation of planar vertex set.

	Pre:	v	List of the interior vertices.

	Post:	Returns vertex data and incidence relations of the

			Ouput: List of sublists {v, f}. Each element of
v is a vertex in
			the triangulation. Each element of f is a list
of three numbers,
			which are the indices of elements in v, defining
a triangle. Thus
			f[[k]] = {a, b, c}, where a, b, c "point" to
			v[[a]], v[[b]], v[[c]].

	Rem:	The adjacency sublists "a" returned by
"DelauneyTriangulation are
			not cyclic. Hence, if a[[1]] is an interior
point, it is necessary
			to add the last face incident a[[1]] before
deleting vertices from
			the triangulation.

Module[ { i, m, n, t = {}, f = {} },

(*	Triangulate vertex list *)

	t = DelauneyTriangulation[v];

(*	Insert cyclic adjacency for t[[i, 1]] if it is interior. *)

	For[i = 1, i <= Length[t], i++,
		If[	MemberQ[ t[[ t[[i, 2, 1]], 2 ]], t[[i, 2, -1]]
			Append[ t[[i, 2]], t[[i, 2, 1]] ]

(*	Build list of faces from t *)

	While[Length[t] > 0,
			Map[Prepend[#, t[[1, 1]]]&, Partition[t[[1, 2]],
2, 1]]

		t = Select[t, ( MemberQ[ Flatten[t[[1]]], #[[1]] ] ==
False ) &];

	{v, f}



We begin with a list of vertices, which might look like this:

    v = {..., v[[a]], ..., v[[p]], ..., v[[q]], ... }
where v[[a]], etc., are points in the plane.

A triangle can be specified by three indices, such as {a, p, q}.

The command "t = DelauneyTriangulation[v]" returns a list of the form:

    t = {{a, {p, q, r, s, t, u}}, ..., {p, {z, q, a, u, x, y}}, ... }

where all entries a, p, q, etc. are "labels" for vertices. The labels
by "DelauneyTriangulation" are the positions of elements in the vertex
(if I remember correctly). I.e., the label for v[[a]] is "a". Let's call
t[[i, 1]] the "anchor" of the adjacency list t[[i]].

The purpose of the "Triangulate" function is to turn the list "t" into
something of the form

    f = {{a, p, q}, {a, q, r}, ..., {a, u, p}, ..., {p, z, q}, {p, q,
a}, ...}

Where each element of f represents a triangle in the Delauney
triangulation t
of v by the labels of it's vertices.

First, it is convenient to augment the adjacency lists in t to make them
"completely cyclic" for interior points.

Consider the following diagram, for example (imagine equilateral

  r   q   z

s   a   p   y
  t   u   x
The vertex v[[a]] will be an interior point of the triangulation if {u,
defines the edge of a triangle, in this case, {a, p, u}.

Description of "For" loop following the comment:

(*	Insert cyclic adjacency for t[[i, 1]] if it is interior. *)

Let i = 1. Then, using our example, above:

t[[i, 1]] = a

t[[ t[[i, 2, 1]], 2 ]] = t[[p, 2]] = {z, q, a, u, x, y}

t[[i, 2, -1]] = u


    MemberQ[ t[[ t[[i, 2, 1]], 2 ]], t[[i, 2, -1]] ]
	    = MemberQ[ {z, q, a, u, x, y}, u ]

which is true. Hence we make the change

	t[[i]] = {a, {p, q, r, s, t, u}} -> {a, {p, q, r, s, t, u, p}}

Now increment "i" and loop.

Description of "While "loop following the comment:

(*	Build list of faces from t *)

First, take t[[1]] = {a, {p, q, r, s, t, u, a}}, and partition it
as follows:

    {p, q, r, s, t, u, a} -> {{p, q}, {q, r}, ... {u, p}}

Now prepend "a" to each sublist to get the list of triangles incident
to the vertex "a", which we will call Star(a):

    Star(a) = {{a, p, q}, {a, q, r}, ... {a, u, p}}

Now join this list to the list of faces "f" (which starts out empty).

Recall that "t" is a list of sublists. The statement

	t = Select[t, ( MemberQ[ Flatten[t[[1]]], #[[1]] ] == False ) &]

selects every element from "t" whose first entry is not in

	{a, p, q, r, s, t, u, p}

Repeat this process until "t" is empty.

The resulting list "f" will have a unique representative for each
face in the triangulation.

To see why this is true, note that each triangle has three
representations in
the augmented adjacency list "t". For example, referring to the above
consider the triangle {p, z, q}, and note that z is not in
{a, p, q, r, s, t, u, p}. The adjacency lists for p, q, z, in "t" are of

    {p, {..., z, ..., q, ...}}
    {q, {..., p, ..., z, ...}}
    {z, {..., q, ..., p, ...}}

The "Select" statement above will drop the first two of these, but
retain the
third, so we still get a representative of triangle {z, p, y} in a later
I claim this always works, for boundary vertices as well as interior

The above process effectively deletes all adjacency lists from t whose
anchor is a vertex of Star(a). Now I will try to justify the assertion
every triangle in the Delauney triangulation is faithfully recorded in
list f when the process is finished. Consider the following cases:

Case 1. The edge {q, p} is interior, and must be incident to two
Hence there is a vertex z, not incident to Star(a), such that {q, p, z}
in Star(z). Since z is not incident to Star(a), the adjacency list with
anchor z survives the deletion process, and the face {q, p, z} is
in a future pass.

Case 2. The edge {q, p} is on the boundary, and each vertex in {q, p} is
incident to at least two edges. In this case, the only triangle incident
{q, p} is {a, p, q} in Star(a).

The adjacency list with anchor p must contain a vertex y, which is not
Star(a). Thus Star(y) contains a triangle which is incident to Star(a)
at vertex p, but not q. Thus, any triangle incident to p, but not in
Star(a) has a representative in t after deletion, and will be recorded
as a face later. (By symmetry, this argument applies to q as well.)

g   h   i

  r   q

s   a   p   y
  t   u   x

Case 3. The edge {q, p} is on the boundary, and one vertex, say p, is
isolated. Then p is only adjacent to a and q, and Star(p) is a single
triangle, which is contained in Star(a). Again, the face {a, p, q} is
recorded once, and the adjacency list with anchor p is deleted without

g   h   i

  r   q

s   a   p

  • Prev by Date: Re: including files
  • Next by Date: Slowdown
  • Previous by thread: DelaunayTriangulation
  • Next by thread: Modeling reactions