Jan Brohl avatar Jan Brohl committed 5d6ff19

sourceformat, trianglefs ->list

Comments (0)

Files changed (1)

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 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 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)
 
-# 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)
 
-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)  
+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 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
-  
- 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()  
+    def flip(self):
+        self.p.reverse()
+
+# return triangle as an ascii STL facet
+    def __str__(self):
+        return ("facet normal 0 0 0\nouter loop\nvertex %s\nvertex %s\nvertex %s\nendloop\nendfacet"
+                % (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)
+
+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)
+
+
+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_dict = {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}
+trianglefs = list(trianglefs_dict[i] for i in range(16))
+
+
+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()
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.