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

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

Python, 302 lines

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.

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, 10 months 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 ,

Alexander Pletzer (author) 8 years, 10 months ago

Besho,

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

--Alex

besho 8 years, 10 months 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

Alexander Pletzer (author) 8 years, 10 months ago

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])

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

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

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

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, 6 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 7 years, 11 months 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) 7 years, 11 months 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 7 years, 7 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) 7 years, 7 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 7 years, 6 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) 7 years, 6 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)