Welcome, guest | Sign In | My Account | Store | Cart

Given a set of points, this code will construct a Delaunay triangulation.

Python, 302 lines
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
import numpy
import math
import copy

class Delaunay2d:

  EPS = 1.23456789e-14

  def __init__(self, points):

    # data structures
    self.points = points[:] # copy
    self.triangles = [] # cells
    self.edge2Triangles = {} # edge to triangle(s) map
    self.boundaryEdges = set()
    self.appliedBoundaryEdges = None
    self.holes = None

    
    # compute center of gravity
    cg = numpy.zeros((2,), numpy.float64)
    for pt in points:
      cg += pt
    cg /= len(points)

    # sort
    def distanceSquare(pt):
      d = pt - cg
      return numpy.dot(d, d)
    self.points.sort(key = distanceSquare)

    # create first triangle, make sure we're getting a non-zero area otherwise
    # drop the points
    area = 0.0
    index = 0
    stop = False
    while not stop and index + 2 < len(points):
      area = self.getArea(index, index + 1, index + 2)
      if abs(area) < self.EPS:
        del self.points[index]
      else:
        stop = True
    if index <= len(self.points) - 3:
      tri = [index, index + 1, index + 2]
      self.makeCounterClockwise(tri)
      self.triangles.append(tri)

      # boundary edges
      e01 = (tri[0], tri[1])
      self.boundaryEdges.add(e01)
      e12 = (tri[1], tri[2])
      self.boundaryEdges.add(e12)
      e20 = (tri[2], tri[0])
      self.boundaryEdges.add(e20)

      e01 = self.makeKey(e01[0], e01[1])
      self.edge2Triangles[e01] = [0,]

      e12 = self.makeKey(e12[0], e12[1])
      self.edge2Triangles[e12] = [0,]

      e20 = self.makeKey(e20[0], e20[1])
      self.edge2Triangles[e20] = [0,]

    else:
      # all the points fall on a line
      return

    # add additional points
    for i in range(3, len(self.points)):
      self.addPoint(i)

    # remove all triangles inside holes
    # TO DO 

  def getTriangles(self):
    """
    @return triangles
    """
    return self.triangles

  def getEdges(self):
    """
    @return egdes
    """
    return self.edge2Triangles.keys()

  def getArea(self, ip0, ip1, ip2):
    """
    Compute the parallelipiped area
    @param ip0 index of first vertex
    @param ip1 index of second vertex
    @param ip2 index of third vertex
    """
    d1 = self.points[ip1] - self.points[ip0]
    d2 = self.points[ip2] - self.points[ip0]
    return (d1[0]*d2[1] - d1[1]*d2[0])

  def isEdgeVisible(self, ip, edge):
    """
    Return true iff the point lies to its right when the edge points down
    @param ip point index
    @param edge (2 point indices with orientation)
    @return True if visible    
    """
    area = self.getArea(ip, edge[0], edge[1])
    if area < self.EPS:
      return True
    return False

  def makeCounterClockwise(self, ips):
    """
    Re-order nodes to ensure positive area (in-place operation)
    """
    area = self.getArea(ips[0], ips[1], ips[2])
    if area < -self.EPS:
      ip1, ip2 = ips[1], ips[2]
      # swap
      ips[1], ips[2] = ip2, ip1

  def flipOneEdge(self, edge):
    """
    Flip one edge then update the data structures
    @return set of edges to interate over at next iteration
    """

    # start with empty set
    res = set()

    # assume edge is sorted
    tris = self.edge2Triangles.get(edge, [])
    if len(tris) < 2:
        # nothing to do, just return
        return res

    iTri1, iTri2 = tris
    tri1 = self.triangles[iTri1]
    tri2 = self.triangles[iTri2]

    # find the opposite vertices, not part of the edge
    iOpposite1 = -1
    iOpposite2 = -1
    for i in range(3):
      if not tri1[i] in edge:
        iOpposite1 = tri1[i]
      if not tri2[i] in edge:
        iOpposite2 = tri2[i]

    # compute the 2 angles at the opposite vertices
    da1 = self.points[edge[0]] - self.points[iOpposite1]
    db1 = self.points[edge[1]] - self.points[iOpposite1]
    da2 = self.points[edge[0]] - self.points[iOpposite2]
    db2 = self.points[edge[1]] - self.points[iOpposite2]
    crossProd1 = self.getArea(iOpposite1, edge[0], edge[1])
    crossProd2 = self.getArea(iOpposite2, edge[1], edge[0])
    dotProd1 = numpy.dot(da1, db1)
    dotProd2 = numpy.dot(da2, db2)
    angle1 = abs(math.atan2(crossProd1, dotProd1))
    angle2 = abs(math.atan2(crossProd2, dotProd2))
    
    # Delaunay's test
    if angle1 + angle2 > math.pi*(1.0 + self.EPS):

      # flip the triangles
      #             / ^ \                    / b \
      # iOpposite1 + a|b + iOpposite2  =>   + - > +
      #             \   /                    \ a /

      newTri1 = [iOpposite1, edge[0], iOpposite2] # triangle a
      newTri2 = [iOpposite1, iOpposite2, edge[1]] # triangle b

      # update the triangle data structure
      self.triangles[iTri1] = newTri1
      self.triangles[iTri2] = newTri2

      # now handle the topolgy of the edges

      # remove this edge
      del self.edge2Triangles[edge]

      # add new edge
      e = self.makeKey(iOpposite1, iOpposite2)
      self.edge2Triangles[e] = [iTri1, iTri2]

      # modify two edge entries which now connect to 
      # a different triangle
      e = self.makeKey(iOpposite1, edge[1])
      v = self.edge2Triangles[e]
      for i in range(len(v)):
        if v[i] == iTri1:
          v[i] = iTri2
      res.add(e)

      e = self.makeKey(iOpposite2, edge[0])
      v = self.edge2Triangles[e]
      for i in range(len(v)):
        if v[i] == iTri2:
          v[i] = iTri1
      res.add(e)

      # these two edges might need to be flipped at the
      # next iteration
      res.add(self.makeKey(iOpposite1, edge[0]))
      res.add(self.makeKey(iOpposite2, edge[1]))

    return res 

  def flipEdges(self):
    """
    Flip edges to statisfy Delaunay's criterion
    """

    # start with all the edges
    edgeSet = set(self.edge2Triangles.keys())

    continueFlipping = True

    while continueFlipping:

      #
      # iterate until there are no more edges to flip
      #

      # collect the edges to flip
      newEdgeSet = set()
      for edge in edgeSet:
        # union
        newEdgeSet |= self.flipOneEdge(edge)

      edgeSet = copy.copy(newEdgeSet)
      continueFlipping = (len(edgeSet) > 0)

  def addPoint(self, ip):
    """
    Add point
    @param ip point index
    """

    # collection for later updates
    boundaryEdgesToRemove = set()
    boundaryEdgesToAdd = set()

    for edge in self.boundaryEdges:

      if self.isEdgeVisible(ip, edge):

        # create new triangle
        newTri = [edge[0], edge[1], ip]
        newTri.sort()
        self.makeCounterClockwise(newTri)
        self.triangles.append(newTri)

        # update the edge to triangle map
        e = list(edge[:])
        e.sort()
        iTri = len(self.triangles) - 1 
        self.edge2Triangles[tuple(e)].append(iTri)

        # add the two boundary edges
        e1 = [ip, edge[0]]
        e1.sort()
        e1 = tuple(e1)
        e2 = [edge[1], ip]
        e2.sort()
        e2 = tuple(e2)
        v1 = self.edge2Triangles.get(e1, [])
        v1.append(iTri)
        v2 = self.edge2Triangles.get(e2, [])
        v2.append(iTri)
        self.edge2Triangles[e1] = v1
        self.edge2Triangles[e2] = v2

        # keep track of the boundary edges to update
        boundaryEdgesToRemove.add(edge)
        boundaryEdgesToAdd.add( (edge[0], ip) )
        boundaryEdgesToAdd.add( (ip, edge[1]) )

    # update the boundary edges
    for bedge in boundaryEdgesToRemove:
      self.boundaryEdges.remove(bedge)
    for bedge in boundaryEdgesToAdd:
      bEdgeSorted = list(bedge)
      bEdgeSorted.sort()
      bEdgeSorted = tuple(bEdgeSorted)
      if len(self.edge2Triangles[bEdgeSorted]) == 1:
        # only add boundary edge if it does not appear
        # twice in different order
        self.boundaryEdges.add(bedge)


    # recursively flip edges
    flipped = True
    while flipped:
      flipped = self.flipEdges()

  def makeKey(self, i1, i2):
    """
    Make a tuple key such at i1 < i2
    """
    if i1 < i2:
      return (i1, i2)
    return (i2, i1)

The algorithm uses the S-hull method by D A Sinclair (http://www.s-hull.org/paper/s_hull.pdf). The method involves ordering the points in increasing distance from the cloud's center of gravity, creating a triangle with the first three points, and adding the remaining points while contructing triangles between the point and boundary edges - only triangles with a definite sign need to be added (the edge must be visible). Finally, an edge flipping step ensures that the triangles are well formed, i.e. the sum of opposite angles to an edge is < 180 degree (the so-called Delaunay criterion). It takes about 30 seconds to triangulate 1000 points, most of the time is spent in two methods (flipOneEdge and getArea).

Further improvements could include adding support for holes and automatically adding Steiner points. Generally speaking, triangles which have an edge on the boundary (convex hull), tend to be squashed; this problem can be alleviated by adding more points on the boundary. Holes can be treaded by removing triangles, edges, and points that are contained inside a given closed path.

15 comments

Alexander Pletzer (author) 9 years, 1 month ago  # | flag

Example of usage:

import random random.seed(1234) xyPoints = [numpy.array([random.random(), random.random()]) for i in range(1000)] delaunay = Delaunay2d(xyPoints) delaunay.show()

besho 8 years, 11 months ago  # | flag

hey , the code dosen't work as it should be . if you give the class more than 9 points it will create the convex hull but will not use all the points .

please help
besho ,

Alexander Pletzer (author) 8 years, 11 months ago  # | flag

Besho,

If you have a driver that exhibits the problem I can take a look at it.

--Alex

besho 8 years, 11 months ago  # | flag

Alex ,

the thing is that when i use your example of usage it dosent work , there is a missing function show() .further more the triangles intersect with each other .

thnx a lot besho

Alexander Pletzer (author) 8 years, 10 months ago  # | flag

The code plot is:

def show(self, width=500, height=500, showVertices=False, showCells=False, showContour=[]):

import Tkinter

xmin = min([p[0] for p in self.points])
ymin = min([p[1] for p in self.points])
xmax = max([p[0] for p in self.points])
ymax = max([p[1] for p in self.points])
padding = 5
w = width - 2*padding
h = height - 2*padding

master = Tkinter.Tk()
c = Tkinter.Canvas(master, width=width, height=height)
c.pack()
for e in self.edge2Triangles:
  i1, i2 = e
  xp1 = padding + int(w*(self.points[i1][0] - xmin)/(xmax - xmin))
  yp1 = padding + int(h*(ymax - self.points[i1][1])/(ymax - ymin))
  xp2 = padding + int(w*(self.points[i2][0] - xmin)/(xmax - xmin))
  yp2 = padding + int(h*(ymax - self.points[i2][1])/(ymax - ymin))
  c.create_line(xp1, yp1, xp2, yp2)

if showVertices:
  for i in range(len(self.points)):
    xp = padding + int(w*(self.points[i][0] - xmin)/(xmax - xmin))
    yp = padding + int(h*(ymax - self.points[i][1])/(ymax - ymin))
    c.create_text(xp, yp, text=str(i))

if showCells:
  for tId, tVals in self.triangles.items():
    cg = reduce(operator.add, [self.points[i] for i in tVals])/float(len(tVals))
    xp = padding + int(w*(cg[0] - xmin)/(xmax - xmin))
    yp = padding + int(h*(ymax - cg[1])/(ymax - ymin))
    c.create_text(xp, yp, text=str(tId))

if len(showContour) > 0:
  for i in range(len(showContour) - 1):
    xp1 = padding + int(w*(showContour[i][0] - xmin)/(xmax - xmin))
    yp1 = padding + int(h*(ymax - showContour[i][1])/(ymax - ymin))
    xp2 = padding + int(w*(showContour[i+1][0] - xmin)/(xmax - xmin))
    yp2 = padding + int(h*(ymax - showContour[i+1][1])/(ymax - ymin))
    c.create_line(xp1, yp1, xp2, yp2, fill='red')


Tkinter.mainloop()
besho 8 years, 10 months ago  # | flag

thnx a lot man its mean alot :) .

if you will ran this code for the plot part below you will see that there is something that dosent tringulate all the points that is given , maybe the version what you have and what is uploaded here is the same .

thnx again man i appreciate it besho ,

import pylab import random

random.seed(50) xyPoints = [numpy.array([random.random(), random.random()]) for i in range(20)] delaunay = Delaunay2d(xyPoints)

xyArray = numpy.asarray(xyPoints)

print xyArray.shape[0], len(delaunay.GetTriangles) , len(delaunay.GetEdges)

for t in delaunay.GetTriangles: # t[0], t[1], t[2] are the points indexes of the triangle t_i = [t[0], t[1], t[2], t[0]] pylab.plot(xyArray[t_i,0],xyArray[t_i,1])

pylab.plot(xyArray[:,0],xyArray[:,1],'o') pylab.show()

besho 8 years, 10 months ago  # | flag

Alex , sorry dude i found my mistake my bad :) it works perfectly.

a lot of thanx besho

Alexander Pletzer (author) 8 years, 10 months ago  # | flag

Besho,

No worries. Be aware that the vertices get reordered internally. One should not assume that the triangle indexing refers to the input xyPoints, rather they refer to internal representation of points (self.points).

--Alex

afshin 8 years, 7 months ago  # | flag

Hi Alexander, Thank you for sharing the code. Is this code applicable for biological membranes. There we do randomly either move one vertice or flip an edge to reach the minimum curvature energy of the membrane.

bludandelo 7 years, 12 months ago  # | flag

Great piece of work! I was wondering whether a new version of the code is available with the following feature:

1) Holes handling 2) Automatic mesh refinement, so that internal points are added in order to keep the average triangles volume constant throughout the domain.

Thanks!

Alexander Pletzer (author) 7 years, 12 months ago  # | flag

Bludandelo: Regarding holes and refinement, you may want to take a look at pytriangle on github: https://github.com/pletzer/pytriangle/blob/master/README.md.

Afshin: Sorry I'm not familiar with biological membranes. You mention two cases: (a) a vertex gets moved and (b) an edge if flipped because some other constraints than Delaunay are applied. The method flipOneEdge could be modified for (b). Regarding (a), I think you can move vertices freely so long each triangle keeps its positive area, which is defined by the ordering of the vertices. In that case the topology would not change, only the vertices. Now it might be a good idea to call flipEdges when you change a vertex as this might require edge to be flipped but that's your call.

Vismaya 7 years, 8 months ago  # | flag

Alex, I have a set of 2d points and I want to apply delaunay triangulation method to these points. The program above doesn't show any plotted outputs though it is running perfectly. Can you help me to get the plotted graph of the 2d points.?

Alexander Pletzer (author) 7 years, 8 months ago  # | flag

Vismaya, The code above (def show...) shows how to plot using Tkinter. You can use matplotlib, see http://matplotlib.org/examples/pylab_examples/triplot_demo.html for instance.

Gene Dial 7 years, 7 months ago  # | flag

In what environment should this code be run? What version Python? When I try it under Spyder 2.3.9 with Python 3.5.2 I get this error: self.points.sort(key = distanceSquare) \n TypeError: 'key' is an invalid keyword argument for this function.

Alexander Pletzer (author) 7 years, 7 months ago  # | flag

Gene, the code was written in python 2.7 but it should run under 3.5, see below. Save the recipe under the file name delaunay.py. Save the following commands

from delaunay import Delaunay2d
import random 
random.seed(1234)
import numpy
xyPoints = [numpy.array([random.random(), random.random()]) for i in range(10)]
delaunay = Delaunay2d(xyPoints)

in file t.py

$ python --version
Python 3.5.0 :: Anaconda 2.4.1 (64-bit)

Run t.py

$ python t.py 
$

I understand that my version is a little behind yours. You can also run the following check:

Python 3.5.0 |Anaconda 2.4.1 (64-bit)| (default, Sep 13 2015, 10:34:39) 
[GCC 4.4.7 20120313 (Red Hat 4.4.7-1)] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> a = [1,2,3,4,5,6]
>>> def dist(x): return (x - 2.3)**2
... 
>>> a.sort(key=dist)
>>> a
[2, 3, 1, 4, 5, 6]

Does that work for you?