# Commits

committed 3a35e5e

(hopefully) flipping triangles to always face the "outside"
added triangle.flip
split main()

# src/tetrahedrons.py

`-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 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 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`
`+  `
`+ def flip(self):`
`+  self.p.reverse()`
`+  `
`+# 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])      `
` `
`-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  `
`+ dr=lambda x,y,z:data[x][y][z]`
`+ xmin=0`
`+ xmax=len(data)-1`
`+ ymin=0`
`+ ymax=len(data[0])-1`
`+ zmin=0`
`+ zmax=len(data[0][0])-1`
`+ #print(data)`
`+ export(dr,isolevel,xmin,xmax,ymin,ymax,zmin,zmax)`
` `
`-# 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`
`+def export(density_function,isolevel,xmin,xmax,ymin,ymax,zmin,zmax):  `
`+ triangles=[]`
`+ p=[None]*8`
`+ val=[None]*8  `
`+ for i in range(xmin,xmax):  `
`+  for j in range(ymin,ymax):  `
`+   for k in range(zmin,zmax):    `
`+    #print(i,j,k)  `
`+    p[0]=Vector(i,j,k)  `
`+    val[0] = density_function(i,j,k ) `
`+    p[1]=Vector(i+1,j,k)  `
`+    val[1] = density_function(i+1,j,k)  `
`+    p[2]=Vector(i+1,j+1,k)  `
`+    val[2] = density_function(i+1,j+1,k)  `
`+    p[3]=Vector(i,j+1,k)  `
`+    val[3] = density_function(i,j+1,k)  `
`+    p[4]=Vector(i,j,k+1)  `
`+    val[4] = density_function(i,j,k+1)  `
`+    p[5]=Vector(i+1,j,k+1)  `
`+    val[5] = density_function(i+1,j,k+1)  `
`+    p[6]=Vector(i+1,j+1,k+1)  `
`+    val[6] = density_function(i+1,j+1,k+1)  `
`+    p[7]=Vector(i,j+1,k+1)  `
`+    val[7] = density_function(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)  `
` `
`-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]`
`+def export_triangles(triangles): # stl format  `
`+ print("solid points")  `
`+ for tri in triangles:  `
`+  print(tri)  `
`+ print("endsolid points")  `
`    `
`-    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()`
`- `
`+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`
`+  `
`+ l= trianglefs[triindex](g, iso, v0, v1, v2, v3)`
`+ if triindex==7 or (triindex!=8 and triindex&8):`
`+     for tri in l:`
`+         tri.flip()`
`+ return l`
`+  `
`+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()  `