Commits

Jan Brohl committed 8b15c22

make it render in panda

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
+from panda3d.core import Vec3, MeshDrawer # @UnresolvedImport
 
-    def __str__(self):
-        return str(self.x) + " " + str(self.y) + " " + str(self.z)
+from math import cos, exp, atan2
+
+
+def lobes(x, y, z):
+    try:
+        theta = atan2(x, y)         # sin t = 0
+    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)
 
 
 class Gridcell:  # struct GRIDCELL
     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]))
+    def __repr__(self):
+        return "%s(%r, %r, %r)" % (self.__class__.__name__,
+                                   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):
+    size = float(size)
     m = int(steps / 2)
     ki = []
     for i in range(steps):
         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 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 Vec3(p1.x + mu * (p2.x - p1.x), p1.y + mu * (p2.y - p1.y), p1.z + mu * (p2.z - p1.z))
 
 
 def t000F(g, iso, v0, v1, v2, v3):
     ]
 
 
-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_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  
+    #   Determine which of the 16 cases we have given which vertices
+    #   are above or below the isosurface
 
     triindex = 0
     if g.val[v0] < iso:
     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))
+def make(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] = Vec3(i, j, k)
+                val[0] = density_function(i, j, k)
+                p[1] = Vec3(i + 1, j, k)
+                val[1] = density_function(i + 1, j, k)
+                p[2] = Vec3(i + 1, j + 1, k)
+                val[2] = density_function(i + 1, j + 1, k)
+                p[3] = Vec3(i, j + 1, k)
+                val[3] = density_function(i, j + 1, k)
+                p[4] = Vec3(i, j, k + 1)
+                val[4] = density_function(i, j, k + 1)
+                p[5] = Vec3(i + 1, j, k + 1)
+                val[5] = density_function(i + 1, j, k + 1)
+                p[6] = Vec3(i + 1, j + 1, k + 1)
+                val[6] = density_function(i + 1, j + 1, k + 1)
+                p[7] = Vec3(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))
+
+    return triangles
+
 
 if __name__ == "__main__":
-    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)
+    import direct.directbase.DirectStart
+    tris = make(dr, isolevel, xmin, xmax, ymin, ymax, zmin, zmax)
+    md = MeshDrawer()
+    mdNode = md.getRoot()
+    mdNode.setTwoSided(True)
+    mdNode.reparentTo(render)  # @UndefinedVariable
+    md.setBudget(xmax * ymax * zmax * 1) #actually 12 instead of 1 in general
+    md.begin(base.cam, render)  # @UndefinedVariable
+    c = (1, 1, 1, 0.5)
+    uv = (0, 0)
+    for tri in tris:
+        v1, v2, v3 = tri.p
+        md.tri(v1, c, uv, v2, c, uv, v3, c, uv)
+    md.end()
+    run()  # @UndefinedVariable