Marching Tetrahedrons in Python

In this atricle we show a simple implementation of the Marching Tetrahedrons algorithm in Python.

Marching Tetrahedrons

The listing povided below is a straightforward reimplementation in Python of the ideas and code presented by Paul Bourke. The image

shows the result of sampling a simple lobed function (rendered in Blender). I might implement the code directly in Blender but for now we export to an STL file that can be read by almost any 3D package. Note that the exported triangles do not have normals that point in a uniform direction, in fact we do no export any normals at all. You have to recalculate the normals in your 3D package before you can render the resulte with smooth shading.

view plainprint?

1. class Vector: # struct XYZ  
2.  def __init__(self,x,y,z):  
3.   self.x=x  
4.   self.y=y  
5.   self.z=z  
6.    
7.  def __str__(self):  
8.   return str(self.x)+" "+str(self.y)+" "+str(self.z)  
9.     
10. class Gridcell: # struct GRIDCELL  
11.  def __init__(self,p,n,val):  
12.   self.p   = p   # p=[8]  
13.   self.n   = n   # n=[8]  
14.   self.val = val # val=[8]  
15.   
16. class Triangle: # struct TRIANGLE  
17.  def __init__(self,p1,p2,p3):  
18.   self.p = [p1, p2, p3] # vertices  
19.   
20. # return triangle as an ascii STL facet    
21.  def __str__(self):  
22.   return """facet normal 0 0 0 
23. outer loop 
24. vertex %s 
25. vertex %s 
26. vertex %s 
27. endloop 
28. endfacet"""%(self.p[0],self.p[1],self.p[2])  
29.   
30. # return a 3d list of values  
31. def readdata(f=lambda x,y,z:x*x+y*y+z*z,size=5.0,steps=11):  
32.  m=int(steps/2)  
33.  ki = []  
34.  for i in range(steps):  
35.   kj = []  
36.   for j in range(steps):  
37.    kd=[]  
38.    for k in range(steps):  
39.     kd.append(f(size*(i-m)/m,size*(j-m)/m,size*(k-m)/m))  
40.    kj.append(kd)  
41.   ki.append(kj)  
42.  return ki  
43.   
44. from math import cos,exp,atan2  
45.   
46. def lobes(x,y,z):  
47.  try:  
48.   theta = atan2(x,y)         # sin t = o   
49.  except:  
50.   theta = 0  
51.  try:  
52.   phi = atan2(z,y)  
53.  except:  
54.   phi = 0  
55.  r = x*x+y*y+z*z  
56.  ct=cos(theta)  
57.  cp=cos(phi)  
58.  return ct*ct*cp*cp*exp(-r/10)  
59.    
60. def main():  
61.   
62.  data = readdata(lobes,5,41)  
63.  isolevel = 0.1  
64.   
65.  #print(data)  
66.    
67.  triangles=[]  
68.  for i in range(len(data)-1):  
69.   for j in range(len(data[i])-1):  
70.    for k in range(len(data[i][j])-1):  
71.     p=[None]*8  
72.     val=[None]*8  
73.     #print(i,j,k)  
74.     p[0]=Vector(i,j,k)  
75.     val[0] = data[i][j][k]  
76.     p[1]=Vector(i+1,j,k)  
77.     val[1] = data[i+1][j][k]  
78.     p[2]=Vector(i+1,j+1,k)  
79.     val[2] = data[i+1][j+1][k]  
80.     p[3]=Vector(i,j+1,k)  
81.     val[3] = data[i][j+1][k]  
82.     p[4]=Vector(i,j,k+1)  
83.     val[4] = data[i][j][k+1]  
84.     p[5]=Vector(i+1,j,k+1)  
85.     val[5] = data[i+1][j][k+1]  
86.     p[6]=Vector(i+1,j+1,k+1)  
87.     val[6] = data[i+1][j+1][k+1]  
88.     p[7]=Vector(i,j+1,k+1)  
89.     val[7] = data[i][j+1][k+1]  
90.      
91.     grid=Gridcell(p,[],val)  
92.     triangles.extend(PolygoniseTri(grid,isolevel,0,2,3,7))  
93.     triangles.extend(PolygoniseTri(grid,isolevel,0,2,6,7))  
94.     triangles.extend(PolygoniseTri(grid,isolevel,0,4,6,7))  
95.     triangles.extend(PolygoniseTri(grid,isolevel,0,6,1,2))  
96.     triangles.extend(PolygoniseTri(grid,isolevel,0,6,1,4))  
97.     triangles.extend(PolygoniseTri(grid,isolevel,5,6,1,4))  
98.       
99.  export_triangles(triangles)  
100.   
101. def export_triangles(triangles): # stl format  
102.  print("solid points")  
103.  for tri in triangles:  
104.   print(tri)  
105.  print("endsolid points")  
106.    
107. def t000F(g, iso, v0, v1, v2, v3):  
108.  return []  
109.   
110. def t0E01(g, iso, v0, v1, v2, v3):  
111.  return [Triangle(  
112.  VertexInterp(iso,g.p[v0],g.p[v1],g.val[v0],g.val[v1]),  
113.  VertexInterp(iso,g.p[v0],g.p[v2],g.val[v0],g.val[v2]),  
114.  VertexInterp(iso,g.p[v0],g.p[v3],g.val[v0],g.val[v3]))  
115.  ]  
116.   
117. def t0D02(g, iso, v0, v1, v2, v3):  
118.  return [Triangle(  
119.  VertexInterp(iso,g.p[v1],g.p[v0],g.val[v1],g.val[v0]),  
120.  VertexInterp(iso,g.p[v1],g.p[v3],g.val[v1],g.val[v3]),  
121.  VertexInterp(iso,g.p[v1],g.p[v2],g.val[v1],g.val[v2]))  
122.  ]  
123.   
124. def t0C03(g, iso, v0, v1, v2, v3):  
125.  tri=Triangle(  
126.  VertexInterp(iso,g.p[v0],g.p[v3],g.val[v0],g.val[v3]),  
127.  VertexInterp(iso,g.p[v0],g.p[v2],g.val[v0],g.val[v2]),  
128.  VertexInterp(iso,g.p[v1],g.p[v3],g.val[v1],g.val[v3]))  
129.  return [tri,Triangle(  
130.  tri.p[2],  
131.  VertexInterp(iso,g.p[v1],g.p[v2],g.val[v1],g.val[v2]),  
132.  tri.p[1])  
133.  ]  
134.   
135. def t0B04(g, iso, v0, v1, v2, v3):  
136.  return [Triangle(  
137.  VertexInterp(iso,g.p[v2],g.p[v0],g.val[v2],g.val[v0]),  
138.  VertexInterp(iso,g.p[v2],g.p[v1],g.val[v2],g.val[v1]),  
139.  VertexInterp(iso,g.p[v2],g.p[v3],g.val[v2],g.val[v3]))  
140.  ]  
141.   
142. def t0A05(g, iso, v0, v1, v2, v3):  
143.  tri = Triangle(  
144.  VertexInterp(iso,g.p[v0],g.p[v1],g.val[v0],g.val[v1]),  
145.  VertexInterp(iso,g.p[v2],g.p[v3],g.val[v2],g.val[v3]),  
146.  VertexInterp(iso,g.p[v0],g.p[v3],g.val[v0],g.val[v3]))  
147.  return [tri,Triangle(  
148.  tri.p[0],  
149.  VertexInterp(iso,g.p[v1],g.p[v2],g.val[v1],g.val[v2]),  
150.  tri.p[1])  
151.  ]  
152.   
153. def t0906(g, iso, v0, v1, v2, v3):  
154.  tri=Triangle(  
155.  VertexInterp(iso,g.p[v0],g.p[v1],g.val[v0],g.val[v1]),  
156.  VertexInterp(iso,g.p[v1],g.p[v3],g.val[v1],g.val[v3]),  
157.  VertexInterp(iso,g.p[v2],g.p[v3],g.val[v2],g.val[v3]))  
158.  return [tri,  
159.  Triangle(  
160.  tri.p[0],  
161.  VertexInterp(iso,g.p[v0],g.p[v2],g.val[v0],g.val[v2]),  
162.  tri.p[2])  
163.  ]  
164.   
165. def t0708(g, iso, v0, v1, v2, v3):  
166.  return [Triangle(  
167.  VertexInterp(iso,g.p[v3],g.p[v0],g.val[v3],g.val[v0]),  
168.  VertexInterp(iso,g.p[v3],g.p[v2],g.val[v3],g.val[v2]),  
169.  VertexInterp(iso,g.p[v3],g.p[v1],g.val[v3],g.val[v1]))  
170.  ]  
171.   
172. trianglefs = {7:t0708,8:t0708,9:t0906,6:t0906,10:t0A05,5:t0A05,11:t0B04,4:t0B04,12:t0C03,3:t0C03,13:t0D02,2:t0D02,14:t0E01,1:t0E01,0:t000F,15:t000F}  
173.   
174. def PolygoniseTri(g, iso, v0, v1, v2, v3):  
175.  triangles = []  
176.   
177.  #   Determine which of the 16 cases we have given which vertices  
178.  #   are above or below the isosurface  
179.   
180.  triindex = 0;  
181.  if g.val[v0] < iso: triindex |= 1  
182.  if g.val[v1] < iso: triindex |= 2  
183.  if g.val[v2] < iso: triindex |= 4  
184.  if g.val[v3] < iso: triindex |= 8  
185.   
186.  return trianglefs[triindex](g, iso, v0, v1, v2, v3)  
187.   
188. def VertexInterp(isolevel,p1,p2,valp1,valp2):  
189.  if abs(isolevel-valp1) < 0.00001 :  
190.   return(p1);  
191.  if abs(isolevel-valp2) < 0.00001 :  
192.   return(p2);  
193.  if abs(valp1-valp2) < 0.00001 :  
194.   return(p1);  
195.  mu = (isolevel - valp1) / (valp2 - valp1)  
196.  return Vector(p1.x + mu * (p2.x - p1.x), p1.y + mu * (p2.y - p1.y), p1.z + mu * (p2.z - p1.z))  
197.   
198. if __name__ == "__main__":  
199.  main()  
200.    

class Vector: # struct XYZ
 def __init__(self,x,y,z):
  self.x=x
  self.y=y
  self.z=z
 
 def __str__(self):
  return str(self.x)+" "+str(self.y)+" "+str(self.z)
  
class Gridcell: # struct GRIDCELL
 def __init__(self,p,n,val):
  self.p   = p   # p=[8]
  self.n   = n   # n=[8]
  self.val = val # val=[8]

class Triangle: # struct TRIANGLE
 def __init__(self,p1,p2,p3):
  self.p = [p1, p2, p3] # vertices

# return triangle as an ascii STL facet  
 def __str__(self):
  return """facet normal 0 0 0
outer loop
vertex %s
vertex %s
vertex %s
endloop
endfacet"""%(self.p[0],self.p[1],self.p[2])

# return a 3d list of values
def readdata(f=lambda x,y,z:x*x+y*y+z*z,size=5.0,steps=11):
 m=int(steps/2)
 ki = []
 for i in range(steps):
  kj = []
  for j in range(steps):
   kd=[]
   for k in range(steps):
    kd.append(f(size*(i-m)/m,size*(j-m)/m,size*(k-m)/m))
   kj.append(kd)
  ki.append(kj)
 return ki

from math import cos,exp,atan2

def lobes(x,y,z):
 try:
  theta = atan2(x,y)         # sin t = o 
 except:
  theta = 0
 try:
  phi = atan2(z,y)
 except:
  phi = 0
 r = x*x+y*y+z*z
 ct=cos(theta)
 cp=cos(phi)
 return ct*ct*cp*cp*exp(-r/10)
 
def main():

 data = readdata(lobes,5,41)
 isolevel = 0.1

 #print(data)
 
 triangles=[]
 for i in range(len(data)-1):
  for j in range(len(data[i])-1):
   for k in range(len(data[i][j])-1):
    p=[None]*8
    val=[None]*8
    #print(i,j,k)
    p[0]=Vector(i,j,k)
    val[0] = data[i][j][k]
    p[1]=Vector(i+1,j,k)
    val[1] = data[i+1][j][k]
    p[2]=Vector(i+1,j+1,k)
    val[2] = data[i+1][j+1][k]
    p[3]=Vector(i,j+1,k)
    val[3] = data[i][j+1][k]
    p[4]=Vector(i,j,k+1)
    val[4] = data[i][j][k+1]
    p[5]=Vector(i+1,j,k+1)
    val[5] = data[i+1][j][k+1]
    p[6]=Vector(i+1,j+1,k+1)
    val[6] = data[i+1][j+1][k+1]
    p[7]=Vector(i,j+1,k+1)
    val[7] = data[i][j+1][k+1]
   
    grid=Gridcell(p,[],val)
    triangles.extend(PolygoniseTri(grid,isolevel,0,2,3,7))
    triangles.extend(PolygoniseTri(grid,isolevel,0,2,6,7))
    triangles.extend(PolygoniseTri(grid,isolevel,0,4,6,7))
    triangles.extend(PolygoniseTri(grid,isolevel,0,6,1,2))
    triangles.extend(PolygoniseTri(grid,isolevel,0,6,1,4))
    triangles.extend(PolygoniseTri(grid,isolevel,5,6,1,4))
    
 export_triangles(triangles)

def export_triangles(triangles): # stl format
 print("solid points")
 for tri in triangles:
  print(tri)
 print("endsolid points")
 
def t000F(g, iso, v0, v1, v2, v3):
 return []

def t0E01(g, iso, v0, v1, v2, v3):
 return [Triangle(
 VertexInterp(iso,g.p[v0],g.p[v1],g.val[v0],g.val[v1]),
 VertexInterp(iso,g.p[v0],g.p[v2],g.val[v0],g.val[v2]),
 VertexInterp(iso,g.p[v0],g.p[v3],g.val[v0],g.val[v3]))
 ]

def t0D02(g, iso, v0, v1, v2, v3):
 return [Triangle(
 VertexInterp(iso,g.p[v1],g.p[v0],g.val[v1],g.val[v0]),
 VertexInterp(iso,g.p[v1],g.p[v3],g.val[v1],g.val[v3]),
 VertexInterp(iso,g.p[v1],g.p[v2],g.val[v1],g.val[v2]))
 ]

def t0C03(g, iso, v0, v1, v2, v3):
 tri=Triangle(
 VertexInterp(iso,g.p[v0],g.p[v3],g.val[v0],g.val[v3]),
 VertexInterp(iso,g.p[v0],g.p[v2],g.val[v0],g.val[v2]),
 VertexInterp(iso,g.p[v1],g.p[v3],g.val[v1],g.val[v3]))
 return [tri,Triangle(
 tri.p[2],
 VertexInterp(iso,g.p[v1],g.p[v2],g.val[v1],g.val[v2]),
 tri.p[1])
 ]

def t0B04(g, iso, v0, v1, v2, v3):
 return [Triangle(
 VertexInterp(iso,g.p[v2],g.p[v0],g.val[v2],g.val[v0]),
 VertexInterp(iso,g.p[v2],g.p[v1],g.val[v2],g.val[v1]),
 VertexInterp(iso,g.p[v2],g.p[v3],g.val[v2],g.val[v3]))
 ]

def t0A05(g, iso, v0, v1, v2, v3):
 tri = Triangle(
 VertexInterp(iso,g.p[v0],g.p[v1],g.val[v0],g.val[v1]),
 VertexInterp(iso,g.p[v2],g.p[v3],g.val[v2],g.val[v3]),
 VertexInterp(iso,g.p[v0],g.p[v3],g.val[v0],g.val[v3]))
 return [tri,Triangle(
 tri.p[0],
 VertexInterp(iso,g.p[v1],g.p[v2],g.val[v1],g.val[v2]),
 tri.p[1])
 ]

def t0906(g, iso, v0, v1, v2, v3):
 tri=Triangle(
 VertexInterp(iso,g.p[v0],g.p[v1],g.val[v0],g.val[v1]),
 VertexInterp(iso,g.p[v1],g.p[v3],g.val[v1],g.val[v3]),
 VertexInterp(iso,g.p[v2],g.p[v3],g.val[v2],g.val[v3]))
 return [tri,
 Triangle(
 tri.p[0],
 VertexInterp(iso,g.p[v0],g.p[v2],g.val[v0],g.val[v2]),
 tri.p[2])
 ]

def t0708(g, iso, v0, v1, v2, v3):
 return [Triangle(
 VertexInterp(iso,g.p[v3],g.p[v0],g.val[v3],g.val[v0]),
 VertexInterp(iso,g.p[v3],g.p[v2],g.val[v3],g.val[v2]),
 VertexInterp(iso,g.p[v3],g.p[v1],g.val[v3],g.val[v1]))
 ]

trianglefs = {7:t0708,8:t0708,9:t0906,6:t0906,10:t0A05,5:t0A05,11:t0B04,4:t0B04,12:t0C03,3:t0C03,13:t0D02,2:t0D02,14:t0E01,1:t0E01,0:t000F,15:t000F}

def PolygoniseTri(g, iso, v0, v1, v2, v3):
 triangles = []

 #   Determine which of the 16 cases we have given which vertices
 #   are above or below the isosurface

 triindex = 0;
 if g.val[v0] < iso: triindex |= 1
 if g.val[v1] < iso: triindex |= 2
 if g.val[v2] < iso: triindex |= 4
 if g.val[v3] < iso: triindex |= 8

 return trianglefs[triindex](g, iso, v0, v1, v2, v3)

def VertexInterp(isolevel,p1,p2,valp1,valp2):
 if abs(isolevel-valp1) < 0.00001 :
  return(p1);
 if abs(isolevel-valp2) < 0.00001 :
  return(p2);
 if abs(valp1-valp2) < 0.00001 :
  return(p1);
 mu = (isolevel - valp1) / (valp2 - valp1)
 return Vector(p1.x + mu * (p2.x - p1.x), p1.y + mu * (p2.y - p1.y), p1.z + mu * (p2.z - p1.z))

if __name__ == "__main__":
 main()