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, tri) self.boundaryEdges.add(e01) e12 = (tri, tri) self.boundaryEdges.add(e12) e20 = (tri, tri) self.boundaryEdges.add(e20) e01 = self.makeKey(e01, e01) self.edge2Triangles[e01] = [0,] e12 = self.makeKey(e12, e12) self.edge2Triangles[e12] = [0,] e20 = self.makeKey(e20, e20) 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*d2 - d1*d2) 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, edge) 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, ips, ips) if area < -self.EPS: ip1, ip2 = ips, ips # swap ips, ips = 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] - self.points[iOpposite1] db1 = self.points[edge] - self.points[iOpposite1] da2 = self.points[edge] - self.points[iOpposite2] db2 = self.points[edge] - self.points[iOpposite2] crossProd1 = self.getArea(iOpposite1, edge, edge) crossProd2 = self.getArea(iOpposite2, edge, edge) 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, iOpposite2] # triangle a newTri2 = [iOpposite1, iOpposite2, edge] # 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) 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) 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)) res.add(self.makeKey(iOpposite2, edge)) 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, edge, 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] e1.sort() e1 = tuple(e1) e2 = [edge, 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, ip) ) boundaryEdgesToAdd.add( (ip, edge) ) # 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. Alexander Pletzer (author) 7 years, 3 months ago

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 7 years ago

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 .

besho , Besho,

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

--Alex besho 7 years ago

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 The code plot is:

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

``````import Tkinter

xmin = min([p for p in self.points])
ymin = min([p for p in self.points])
xmax = max([p for p in self.points])
ymax = max([p for p in self.points])

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] - xmin)/(xmax - xmin))
yp1 = padding + int(h*(ymax - self.points[i1])/(ymax - ymin))
xp2 = padding + int(w*(self.points[i2] - xmin)/(xmax - xmin))
yp2 = padding + int(h*(ymax - self.points[i2])/(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] - xmin)/(xmax - xmin))
yp = padding + int(h*(ymax - self.points[i])/(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 - xmin)/(xmax - xmin))
yp = padding + int(h*(ymax - cg)/(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] - xmin)/(xmax - xmin))
yp1 = padding + int(h*(ymax - showContour[i])/(ymax - ymin))
xp2 = padding + int(w*(showContour[i+1] - xmin)/(xmax - xmin))
yp2 = padding + int(h*(ymax - showContour[i+1])/(ymax - ymin))
c.create_line(xp1, yp1, xp2, yp2, fill='red')

Tkinter.mainloop()
`````` besho 7 years ago

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, len(delaunay.GetTriangles) , len(delaunay.GetEdges)

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

pylab.plot(xyArray[:,0],xyArray[:,1],'o') pylab.show() besho 7 years ago

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

a lot of thanx besho 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 6 years, 9 months ago

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 6 years, 1 month ago

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) 6 years, 1 month ago

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 5 years, 10 months ago

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) 5 years, 10 months ago

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 5 years, 9 months ago

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) 5 years, 9 months ago

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
>>> 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? Created by Alexander Pletzer on Sun, 8 Feb 2015 (MIT)