"Simple triangulation algorithm in python"
Posted on 2016-03-17 in programming
Last week for whatever reason I've been tinkering with triangulation algorithms. Basically, I've been attempting to make something similar to marvelous designer's clothing triangulation algorithm.
You see, blender has clothing sim, that is, although quite slow, can be used as a basic replacement for Marvelous Designer. The tools are inferior, of course, but you can make pieces of cloth, put them in midair, and connect to each other with sewing lines so gravity would pull them together. Then you would run a cloth sim, and waste next 30 minutes making sure that cloth sim doesns't explode.
But I digress.
For this kind of tool, you'll need ability to tesselate flat shape into sorta uniform-sized triangles. That is called "constrained delaunay triangulation", and there'ss already triangle library that implement it. Said library could be probably studied and converted to python or whatever other language.
Either way I wanted to play around with tesselation myself first.
The simple algorithm I tried was "find longest edge, subdivide it in half, repeat forever, until there are no edges longer than minimum size".
Here's unsubdivised mesh:
And here's subdivided result:
It doesn't have the uniform structure of marvelous designer's algorithm, but it is a decent start.
And here is the code:
The program is written in pyhon, and requires PyOpenGL to work.
from OpenGL.GLUT import *
from OpenGL.GLU import *
from OpenGL.GL import *
import sys
scrWidth = 800
scrHeight = 600
def findIndex(el, arr):
for i, item in enumerate(arr):
if item == el:
return i
return -1;
class Mesh:
def squareDist(self, idxA, idxB):
vertA = self.verts[idxA]
vertB = self.verts[idxB]
diffX = vertA[0] - vertB[0]
diffY = vertA[1] - vertB[1]
return diffX*diffX + diffY*diffY
def distance(self, idxA, idxB):
return math.sqrt(self.squareDist(idxA, idxB))
def addVertex(self, x, y):
newVert = (x, y)
result = findIndex(newVert, self.verts)
if result >= 0:
return result
result = len(self.verts)
self.verts.append(newVert)
return result
def addFace(self, idxA, idxB, idxC):
newFace = (idxA, idxB, idxC)
faceIndex = self.faceIndex
self.faceIndex += 1
self.faces[faceIndex] = newFace
self.addEdge(idxA, idxB, faceIndex, 0)
self.addEdge(idxB, idxC, faceIndex, 1)
self.addEdge(idxC, idxA, faceIndex, 2)
return faceIndex
def removeFace(self, faceIndex):
foundFace = self.faces.get(faceIndex, None)
if foundFace is None:
raise ValueError("cannot remove face {0}: face does not exist".format(faceIndex))
self.edges.pop((foundFace[0], foundFace[1]))
self.edges.pop((foundFace[1], foundFace[2]))
self.edges.pop((foundFace[2], foundFace[0]))
self.faces.pop(faceIndex, None)
def flipEdge(self, idxA, idxB):
forwardEdge = (idxA, idxB)
reverseEdge = (idxB, idxA)
forwardFaceData = self.edges.get(forwardEdge, None)
reverseFaceData = self.edges.get(reverseEdge, None)
if (forwardFaceData is None) or (reverseFaceData is None):
raise ValueError("Cannot flip edge {0} {1}: not enough faces", idxA, idxB)
forwardFace = self.faces[forwardFaceData[0]]
reverseFace = self.faces[reverseFaceData[0]]
self.removeFace(forwardFaceData[0])
self.removeFace(reverseFaceData[0])
a = forwardFace[forwardFaceData[1]]
b = reverseFace[(reverseFaceData[1] + 2)%3]
c = forwardFace[(forwardFaceData[1] + 1)%3]
d = forwardFace[(forwardFaceData[1] + 2)%3]
self.addFace(d, a, b)
self.addFace(d, b, c)
resultEdge = (d, b)
return resultEdge
def addEdge(self, idxA, idxB, faceIndex, vertexIndex):
key = (idxA, idxB)
foundEdge = self.edges.get(key, None)
if not (foundEdge is None):
raise ValueError("topology error - duplicate edge ({0}, {1}) found".format(idxA, idxB))
self.edges[key] = (faceIndex, vertexIndex, self.squareDist(idxA, idxB))
def __init__(self):
self.faceIndex = 0
self.verts = [] #tuples, (x y) format
self.faces = {} #faces (vertIdx1, vertIdx2, vertIdx3) format
self.edges = {} #edge map (startVert, endVert) -> (faceIndex, faceVertexIndex, length)
pass
def getLongestEdge(self):
result = None
lastDist = 0.0
for key, value in self.edges.iteritems():
newDist = value[2]
if newDist > lastDist:
result = key
lastDist = newDist
return result, lastDist
def subdivideFaceEdge(self, faceIdx, faceEdgeIndex, newVertexIdx):
face = self.faces[faceIdx]
self.removeFace(faceIdx)
idxPrev = face[faceEdgeIndex]
idxNext = face[(faceEdgeIndex + 1)%3]
idxMid = face[(faceEdgeIndex + 2)%3]
self.addFace(idxPrev, newVertexIdx, idxMid)
self.addFace(newVertexIdx, idxNext, idxMid)
pass
def subdivideEdge(self, idxA, idxB):
forwardEdge = (idxA, idxB)
reverseEdge = (idxB, idxA)
forwardFaceData = self.edges.get(forwardEdge, None)
reverseFaceData = self.edges.get(reverseEdge, None)
if (forwardFaceData is None) and (reverseFaceData is None):
raise ValueError("Cannot subdivide: edge not found")
vertA = self.verts[idxA]
vertB = self.verts[idxB]
newVertIndex = self.addVertex(vertA[0]*0.5 + vertB[0]*0.5, vertA[1]*0.5+vertB[1]*0.5)
if not (forwardFaceData is None):
self.subdivideFaceEdge(forwardFaceData[0], forwardFaceData[1], newVertIndex)
if not (reverseFaceData is None):
self.subdivideFaceEdge(reverseFaceData[0], reverseFaceData[1], newVertIndex)
return newVertIndex
def drawVertex(self, idx):
glVertex2f(self.verts[idx][0], self.verts[idx][1])
def drawEdge(self, idx1, idx2):
self.drawVertex(idx1)
self.drawVertex(idx2)
def draw(self):
glBegin(GL_LINES)
for idx, face in self.faces.iteritems():
self.drawEdge(face[0], face[1])
self.drawEdge(face[1], face[2])
self.drawEdge(face[2], face[0])
glEnd()
def main():
glutInit(sys.argv)
glutInitDisplayMode(GLUT_DOUBLE|GLUT_RGB|GLUT_DEPTH)
glutInitWindowSize(800, 600)
glutCreateWindow("subdiv test")
glClearColor(0.,0.,0.,1.)
glShadeModel(GL_SMOOTH)
glEnable(GL_CULL_FACE)
#glEnable(GL_LIGHTING)
glutDisplayFunc(display)
glutKeyboardFunc(keyFunc)
glutMainLoop()
def createMesh():
mesh = Mesh()
a = mesh.addVertex(100., 100.)
b = mesh.addVertex(400., 100.)
c = mesh.addVertex(400., 400.)
d = mesh.addVertex(100., 400.)
e = mesh.addVertex(300., 600.)
f = mesh.addVertex(700., 200.)
g = mesh.addVertex(800., 600.)
mesh.addFace(a, b, c)
mesh.addFace(a, c, d)
mesh.addFace(d, c, e)
mesh.addFace(b, f, c)
mesh.addFace(f, g, c)
mesh.addFace(c, g, e)
#mesh.removeFace(faceIdx)
#newVert = mesh.subdivideEdge(a, c)
#newVert = mesh.subdivideEdge(a, newVert)
#newVert = mesh.subdivideEdge(a, newVert)
#newVert = mesh.subdivideEdge(a, newVert)
#mesh.flipEdge(a, c)
return mesh
def keyFunc(key, x, y):
#print key
if key == " ":
longestEdge = mesh.getLongestEdge()
#print longestEdge
if not (longestEdge[0] is None):
mesh.subdivideEdge(longestEdge[0][0], longestEdge[0][1])
glutPostRedisplay()
if key == "a":
minDist = 25.0
minDist *= minDist
while True:
longestEdge = mesh.getLongestEdge()
print longestEdge
if longestEdge[1] < minDist:
break
if longestEdge[0] is None:
break
mesh.subdivideEdge(longestEdge[0][0], longestEdge[0][1])
glutPostRedisplay()
pass
mesh = createMesh()
def display():
glClear(GL_COLOR_BUFFER_BIT|GL_DEPTH_BUFFER_BIT)
glMatrixMode(GL_PROJECTION)
glLoadIdentity()
gluOrtho2D(0, scrWidth, scrHeight, 0)
glMatrixMode(GL_MODELVIEW)
glLoadIdentity()
mesh.draw()
glutSwapBuffers()
if __name__ == '__main__':main()