RE: DelaunayTriangulation

*To*: mathgroup at smc.vnet.net*Subject*: [mg53210] RE: [mg53081] DelaunayTriangulation*From*: "Melko, Mike" <Mike.Melko at northern.edu>*Date*: Thu, 30 Dec 2004 01:38:15 -0500 (EST)*Sender*: owner-wri-mathgroup at wolfram.com

-----Original Message----- From: peteraptaker [mailto:psa at laplacian.co.uk] To: mathgroup at smc.vnet.net 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 assume). Please excuse by finite element and Dealauney ignorance: Peter ********************************************* Mike Melko wrote: Peter, I'm writing a package that uses "DelauneyTriangulation" for certain purposes, and have spent some time thinking about this. Here is my attempt to answer your questions. Let me say that I have found this a bit confusing, and difficult 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 geometry myself. 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 desired triangulation and boundary nodes (no longer a convex hull I assume). Answer: I actually forgot about this point, and I thank you for calling it to my attention. I believe the answer is that you must be able to identify 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" points, and then deleting those points and their adjacencies from the output of "DelauneyTriangulation". ************************************************************************ ***** 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 data 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 functions to do what you want to do. I wrote a function that does something like what 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 moment) 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 start, anyway. ************************************************************************ ***** (* 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 Delauney triangulation 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 vertices 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. *) Triangulate[v_List]:= 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, Join[ f, Map[Prepend[#, t[[1, 1]]]&, Partition[t[[1, 2]], 2, 1]] ]; t = Select[t, ( MemberQ[ Flatten[t[[1]]], #[[1]] ] == False ) &]; ]; {v, f} ] ************************************************************************ ****** DISESCTION OF FUNCTION 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 used by "DelauneyTriangulation" are the positions of elements in the vertex list (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 triangles): r q z s a p y t u x The vertex v[[a]] will be an interior point of the triangulation if {u, p} 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 Then 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 diagram, 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 the form: {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 pass. I claim this always works, for boundary vertices as well as interior vertices. 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 that every triangle in the Delauney triangulation is faithfully recorded in the 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 triangles. Hence there is a vertex z, not incident to Star(a), such that {q, p, z} is 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 recorded 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 to {q, p} is {a, p, q} in Star(a). The adjacency list with anchor p must contain a vertex y, which is not in 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 harm. g h i r q s a p t