Welcome, guest | Sign In | My Account | Store | Cart
```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])
e12 = (tri[1], tri[2])
e20 = (tri[2], tri[0])

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

for i in range(3, len(self.points)):

# 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
"""

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]

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

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

# these two edges might need to be flipped at the
# next iteration

return res

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

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)

"""
@param ip point index
"""

boundaryEdgesToRemove = 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

# update the boundary edges
for bedge in boundaryEdgesToRemove:
self.boundaryEdges.remove(bedge)
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

# 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)
```