In [1]:
__author__ = "Víctor Adolfo Gallego Alcalá"
import numpy as np
import mpl_toolkits.mplot3d as a3
import matplotlib.colors as colors
import pylab as pl
import scipy as sp
import time
# Classes for the data structures
class Vertex:
""" This class represents a vertex, i.e., a point in the space
Attributes
----------
v : array of float
Coordinates of the points
next : Vertex
Next Vertex in the structure
prev : Vertex
Previous Vertex in the structure
duplicate : Edge
Pointer to incident cone edge
onhull : boolean
True iff vertex is on the hull
mark : boolean
True iff point already processed
name : string
Contains its three coordinates
"""
def __init__(self, x=None, y=None, z=None):
self.v = [x,y,z]
self.next = None
self.prev = None
self.duplicate = None # pointer to incident cone edge
self.onhull = False # True iff point on hull
self.mark = False # True iff point already processed
self.name = "( " + str(x) + ", " + str(y) + ", " + str(z) + " )"
def __repr__(self):
""" Returns a string containing its three coordinates
"""
return self.name
def tolist(self):
""" Returns a list containing its three coordinates
"""
return [self.v[0], self.v[1], self.v[2]]
def tolist2D(self):
""" Returns a list containing its two first coordinates (projection
to z=0)
"""
return [self.v[0], self.v[1]]
class Edge:
""" This class represents an edge, i.e., a segment
connecting two vertices
Attributes
----------
endPts : array of Vertex
The vertices of the segment
faces : array of Face
The two incident faces
prev : Edge
Previous Edge in the structure
next : Edge
Next Edge in the structure
newface : Face
Incident cone face generated
delete : boolean
True iff this object should be deleted
name : string
Contains its two vertices
"""
def __init__(self, v1=None, v2=None):
self.endPts = [v1, v2]
self.faces = [None, None]
self.prev = None
self.next = None
self.newface = None # pointer to incident cone face
self.delete = False # True iff this object should be deleted
if v1 != None and v2 != None:
self.name = v1.name + " - " + v2.name
else:
self.name = "empty edge"
#edges = add(self, edges)
def __repr__(self):
""" Returns a string containing its two vertices
"""
return self.name
class Face:
""" This class represents a (triangular) face
connecting two vertices
Attributes
----------
vertex : array of Vertex
The three vertices of the face
edges : array of Edge
The three edges of the face
prev : Face
Previous Face in the structure
next : Face
Next Face in the structure
newface : Face
Incident cone face generated
visible : boolean
True iff this face is visible from a fixed point
name : string
Contains its three vertices
"""
def __init__(self, v0=None, v1=None, v2=None, foldFace=None, edges=None):
if foldFace == None:
e0, e1, e2 = Edge(), Edge(), Edge()
else:
e0, e1, e2 = reversed(foldFace.edges)
e0.endPts = [v0, v1]
e1.endPts = [v1, v2]
e2.endPts = [v2, v0]
e0.name = "( " + str(v0) + " - " + str(v1) + " )"
e1.name = "( " + str(v1) + " - " + str(v2) + " )"
e2.name = "( " + str(v2) + " - " + str(v0) + " )"
self.edges = [e0, e1, e2]
self.vertex = [v0, v1, v2]
#self.edges = [None, None, None]
self.prev = None
self.next = None
self.visible = False # True iff face visible
# Linking edges to faces
e0.faces[0] = e1.faces[0] = e2.faces[0] = self
if v0 != None and v1 != None and v2 != None:
self.name = v0.name + " - " + v1.name + " - " + v2.name
else:
self.name = None
#add(self, faces)
def __repr__(self):
""" Returns a string containing its three vertices
"""
return self.name
def tolist(self):
""" Returns a list containing its three vertices
"""
return [self.vertex[0].tolist(), self.vertex[1].tolist(), self.vertex[2].tolist()]
def tolist2D(self):
""" Returns a list containing its three vertices projected to z=0
"""
return [self.vertex[0].tolist2D(), self.vertex[1].tolist2D(), self.vertex[2].tolist2D()]
In [2]:
class Context:
""" This class performs the beneath-beyond algorithm
to compute the convex hull of a set of points
in R^3
Attributes
----------
vertices : circular doubly linked list of Vertex
edges : circular doubly linked list of Edge
faces : circular doubly linked list of Face
"""
def __init__(self):
self.vertices = None
self.edges = None
self.faces = None
def addVertex(self, p):
""" Adds a new vertex p to the back of the vertices list
"""
if self.vertices is not None:
p.next = self.vertices
p.prev = self.vertices.prev
self.vertices.prev = p
p.prev.next = p
else:
self.vertices = p
self.vertices.next = self.vertices.prev = p
def addEdge(self, p):
""" Adds a new edge p to the back of the edges list
"""
if self.edges is not None:
p.next = self.edges
p.prev = self.edges.prev
self.edges.prev = p
p.prev.next = p
else:
self.edges = p
self.edges.next = self.edges.prev = p
def addFace(self, p):
""" Adds a new face p to the back of the faces list
"""
if self.faces is not None:
p.next = self.faces
p.prev = self.faces.prev
self.faces.prev = p
p.prev.next = p
else:
self.faces = p
self.faces.next = self.faces.prev = p
def deleteVertex(self, p):
""" Deletes a vertex p from the vertices list
"""
if self.vertices is not None:
if self.vertices is self.vertices.next:
self.vertices = None
elif p is self.vertices:
self.vertices = self.vertices.next
p.next.prev = p.prev
p.prev.next = p.next
def deleteEdge(self, p):
""" Deletes an edge p from the edges list
"""
if self.edges is not None:
if self.edges is self.edges.next:
self.edges = None
elif p is self.edges:
self.edges = self.edges.next
p.next.prev = p.prev
p.prev.next = p.next
def deleteFace(self, p):
""" Deletes a face p from the faces list
"""
if self.faces is not None:
if self.faces is self.faces.next:
self.faces = None
elif p is self.faces:
self.faces = self.faces.next
p.next.prev = p.prev
p.prev.next = p.next
def createVertices(self, L):
""" Given a list of [x,y,z] points it creates the
vertices list
"""
for point in L:
v = Vertex(point[0], point[1], point[2])
self.addVertex(v)
def createVerticesFrom2D(self, L):
""" Similar to createVertices, but with bidimensional
points. The z coordinate is asigned to x^2 + y^2. Used
to compute a Delaunay triangulation
"""
for point in L:
v = Vertex(point[0], point[1], point[0]**2 + point[1]**2)
self.addVertex(v)
def collinear(self, a, b, c):
""" Checks if the three given vertices are collinear, using the cross product
"""
return (c.v[2] - a.v[2])*(b.v[1] - a.v[1]) - (b.v[2] - a.v[2])*(c.v[1] - a.v[1]) == 0 \
and (b.v[2] - a.v[2])*(c.v[0] - a.v[0]) - (b.v[0] - a.v[0])*(c.v[2] - a.v[2]) == 0 \
and (b.v[0] - a.v[0])*(c.v[1] - a.v[1]) - (b.v[1] - a.v[1])*(c.v[0] - a.v[0]) == 0
def doubleTriangle(self):
""" Constructs the initial convex polytope, which consists in the first three
non collinear points and the two triangles they constitute (different orientations
give different triangles). In addition, a fourth point which is non coplanar is found
and set to the head of the vertices list
"""
v0 = self.vertices
while self.collinear(v0, v0.next, v0.next.next):
if v0.next.next is self.vertices:
print "Collinearmania!"
return
v1, v2 = v0.next, v0.next.next
v0.mark = True
v1.mark = True
v2.mark = True
f0 = Face(v0, v1, v2)
f1 = Face(v2, v1, v0, f0)
self.addFace(f0)
self.addFace(f1)
for e in f0.edges:
self.addEdge(e)
# Linking
f0.edges[0].faces[1] = f1
f0.edges[1].faces[1] = f1
f0.edges[2].faces[1] = f1
f1.edges[0].faces[1] = f0
f1.edges[1].faces[1] = f0
f1.edges[2].faces[1] = f0
# Find a fourth, non coplanar point
v3 = v2.next
while (not self.volumeSign(f0, v3)):
if v3 == None:
print "Coplanarmania!"
return
v3 = v3.next
self.vertices = v3
def volumeSign(self, f, p):
""" Computes the sign of the volume of a tetrahedron
formed by the vertices of f and a vertex p. The volume
is computed as a 4x4 determinant
"""
ax = f.vertex[0].v[0] - p.v[0]
ay = f.vertex[0].v[1] - p.v[1]
az = f.vertex[0].v[2] - p.v[2]
bx = f.vertex[1].v[0] - p.v[0]
by = f.vertex[1].v[1] - p.v[1]
bz = f.vertex[1].v[2] - p.v[2]
cx = f.vertex[2].v[0] - p.v[0]
cy = f.vertex[2].v[1] - p.v[1]
cz = f.vertex[2].v[2] - p.v[2]
# Beware of numerical imprecissions!
return np.sign(ax*(by*cz-bz*cy) + ay*(bz*cx - bx*cz) + az*(bx*cy - by*cx))
def constructHull(self):
""" The main function of the algorithm. It adds a vertex
to the convex hull in each iteration, and deletes the hidden vertices,
edges and faces.
"""
v = self.vertices
vnext = v.next
if not v.mark:
v.mark = True
self.addOneVertex(v)
self.cleanUp()
v = vnext
while v is not self.vertices:
vnext = v.next
if not v.mark:
v.mark = True
self.addOneVertex(v)
self.cleanUp()
v = vnext
def addOneVertex(self, p):
""" Adds a point p to the hull, creating a new cone of
faces if p is exterior
"""
vis = False
f = self.faces
# We mark all visible faces
if self.volumeSign(f, p) < 0:
f.visible = True
vis = True
f = f.next
while f is not self.faces:
if self.volumeSign(f, p) < 0:
f.visible = True
vis = True
f = f.next
# If no faces are visibles from p, p is in the hull
if not vis:
p.onhull = False
return False
# Mark edges interior to visible region for future deletion
e = self.edges
tmp = e.next
if e.faces[0].visible and e.faces[1].visible:
e.delete = True
elif e.faces[0].visible or e.faces[1].visible:
# e is a visible border: erect a new face
e.newface = self.makeConeFace(e, p)
e = tmp
while e is not self.edges:
tmp = e.next
if e.faces[0].visible and e.faces[1].visible:
e.delete = True
elif e.faces[0].visible or e.faces[1].visible:
# e is a visible border: erect a new face
e.newface = self.makeConeFace(e, p)
e = tmp
return True
def makeConeFace(self, e, p):
""" Helper method of addOneVertex. It creates a new face
connecting a visible border edge and point p
"""
newEdge = [None]*2
for i in range(2):
newEdge[i] = e.endPts[i].duplicate
if newEdge[i] == None:
newEdge[i] = Edge(e.endPts[i], p)
self.addEdge(newEdge[i])
e.endPts[i].duplicate = newEdge[i]
newFace = Face()
newFace.edges[0] = e
newFace.edges[1] = newEdge[0]
newFace.edges[2] = newEdge[1]
self.addFace(newFace)
self.makeCcw(newFace, e, p)
newFace.name = newFace.vertex[0].name + " - " + newFace.vertex[1].name + " - " + newFace.vertex[2].name
# Linking
for i in range(2):
for j in range(2):
if newEdge[i].faces[j] == None:
newEdge[i].faces[j] = newFace
break
return newFace
def makeCcw(self, f, e, p):
""" Helper method of makeConeFace. It enforces a counter-clockwise
order of the vertices in the new face
"""
if e.faces[0].visible:
fv = e.faces[0]
else:
fv = e.faces[1]
i = 0
while fv.vertex[i] != e.endPts[0]:
i += 1
# Orient f the same as fv
if fv.vertex[ (i+1) % 3 ] != e.endPts[1]:
f.vertex[0] = e.endPts[1]
f.vertex[1] = e.endPts[0]
else:
f.vertex[0] = e.endPts[0]
f.vertex[1] = e.endPts[1]
self.swap(f.edges[1], f.edges[2])
f.vertex[2] = p
def swap(self, A, B):
""" Helper method. It swaps the value
of two variables
"""
tmp = A
A = B
B = tmp
def cleanUp(self):
""" It performs the deletion of old objects
during an iteration of the hull construction
"""
self.cleanEdges()
self.cleanFaces()
self.cleanVertices()
def cleanFaces(self):
""" Removes the non visible faces
"""
while self.faces is not None and self.faces.visible:
f = self.faces
self.deleteFace(f)
f = self.faces.next
if f.visible:
t = f
f = f.next
self.deleteFace(t)
else:
f = f.next
while f is not self.faces:
if f.visible:
t = f
f = f.next
self.deleteFace(t)
else:
f = f.next
def cleanEdges(self):
""" Removes the hidden edges
"""
# Integrate the newfaces into the edges
e = self.edges
if e.newface is not None:
if e.faces[0].visible:
e.faces[0] = e.newface
else:
e.faces[1] = e.newface
e.newface = None
e = e.next
while e is not self.edges:
if e.newface != None:
if e.faces[0].visible:
e.faces[0] = e.newface
else:
e.faces[1] = e.newface
e.newface = None
e = e.next
# Delete marked edges
while self.edges is not None and self.edges.delete:
e = self.edges
self.deleteEdge(e)
e = self.edges.next
if e.delete:
t = e
e = e.next
self.deleteEdge(t)
else:
e = e.next
while e is not self.edges:
if e.delete:
t = e
e = e.next
self.deleteEdge(t)
else:
e = e.next
def cleanVertices(self):
""" Removes the hidden vertices
"""
e = self.edges
e.endPts[0].onhull, e.endPts[1].onhull = True, True
e = e.next
while e is not self.edges:
e.endPts[0].onhull, e.endPts[1].onhull = True, True
e = e.next
while self.vertices is not None and self.vertices.mark and not self.vertices.onhull:
v = self.vertices
self.deleteVertex(v)
v = self.vertices.next
if v.mark and not v.onhull:
t = v
v = v.next
self.deleteVertex(t)
else:
v = v.next
while v is not self.vertices:
if v.mark and not v.onhull:
t = v
v = v.next
self.deleteVertex(t)
else:
v = v.next
v = self.vertices
v.duplicate = None
v.onhull = False
v = v.next
while v is not self.vertices:
v.duplicate = None
v.onhull = True
v = v.next
def normZ(self, f):
""" Computes the z component of the normal to
face f
"""
a = f.vertex[0]
b = f.vertex[1]
c = f.vertex[2]
return (b.v[0] - a.v[0])*(c.v[1] - a.v[1]) - \
(b.v[1] - a.v[1])*(c.v[0] - a.v[0])
def lowerFaces(self):
""" Returns a list of the lower faces in the convex hull, i.e.,
those whose normal points downwards
"""
L = []
f = self.faces
while True:
if self.normZ(f) < 0:
L.append(f)
f = f.next
if f is self.faces:
break
return L
In [4]:
def pr(L):
""" Helper function to print the elemens from any
circular doubly linked list
"""
print "----------------"
print L
v = L.next
n = 0
while v is not L and n < 20:
print (n+1,v)
v = v.next
n += 1
print "----------------"
def plot(tris, lim, data):
""" Plots the 3D convex hull
"""
ax = a3.Axes3D(pl.figure())
ax.set_zlim3d([0, lim])
ax.set_xlim3d([0, lim])
ax.set_ylim3d([0, lim])
ax.scatter3D(data[:,0], data[:,1], data[:,2])
tri = a3.art3d.Poly3DCollection([np.asarray(tris.tolist())])
tri.set_color(colors.rgb2hex(sp.rand(3)))
tri.set_edgecolor('k')
ax.add_collection3d(tri)
v = tris.next
while v is not tris:
tri = a3.art3d.Poly3DCollection([np.asarray(v.tolist())])
tri.set_color(colors.rgb2hex(sp.rand(3)))
tri.set_edgecolor('k')
ax.add_collection3d(tri)
v = v.next
pl.show()
In [10]:
%matplotlib inline
c = Context()
data = sp.rand(70,3)
c.createVertices(data.tolist())
c.doubleTriangle()
c.constructHull()
plot(c.faces, 1, data)
In [11]:
# Delaunay triangulation
import matplotlib.pyplot as plt
c = Context()
data = sp.randn(70,2)
plt.scatter(data[:,0], data[:,1])
c.createVerticesFrom2D(data.tolist())
c.doubleTriangle()
c.constructHull()
L = c.lowerFaces()
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon
plt.scatter(data[:,0], data[:,1])
for i in L:
pts = np.asarray(i.tolist2D())
pts.reshape(3,2,order='C')
p = Polygon(pts, closed=True)
p.set_facecolor(colors.rgb2hex(sp.rand(3)))
p.set_edgecolor('k')
p.set_alpha(1)
ax = plt.gca()
ax.add_patch(p)
plt.show()
In [ ]: