Commits

Jan Brohl  committed 1190039

seperated uncached aproach, COMMENTS

  • Participants
  • Parent commits c1c4d2b

Comments (0)

Files changed (1)

File src/tetrahedrons.py

 from panda3d.core import Vec3, MeshDrawer  # @UnresolvedImport
 from density import Metaball, Composite
+SMALL_DOUBLE = 0.00001
 
 
-class Gridcell(object):  # 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(object):
+    """
+    One cell on the grid
+    """
+    def __init__(self, p, val):
+        self.p = p  # p=[8] points
+        self.val = val  # val=[8] densities
 
 
-class Triangle(object):  # struct TRIANGLE
+class Triangle(object):
+    """
+    A simple Triangle
+    """
     def __init__(self, p1, p2, p3):
         self.p = [p1, p2, p3]  # vertices
 
                                    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):
+def cache(f, size, steps):
+    """
+    Evaluate (density-function) f on a grid and put the results to a list of lists of lists
+    """
     size = float(size)
     m = int(steps / 2)
     ki = []
     return ki
 
 
-def VertexInterp(isolevel, p1, p2, valp1, valp2):
-    if abs(isolevel - valp1) < 0.00001:
+def vertexInterp(isolevel, p1, p2, valp1, valp2):
+    """
+    Linearly interpolate the position where an isosurface cuts
+    an edge between two vertices, each with their own scalar value
+    """
+    if abs(isolevel - valp1) < SMALL_DOUBLE:
         return(p1)
-    if abs(isolevel - valp2) < 0.00001:
+    if abs(isolevel - valp2) < SMALL_DOUBLE:
         return(p2)
-    if abs(valp1 - valp2) < 0.00001:
+    if abs(valp1 - valp2) < SMALL_DOUBLE:
         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))
 
 
+# for t???? see http://paulbourke.net/geometry/polygonise/polytetra2.gif
+# ?? ?? is hex for the two cases in binary notation
+
+
 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]))
+    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]))
+    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]))
+    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]),
+    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]))
+    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]))
+    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]),
+    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]))
+    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]),
+    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]))
+    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]))
     ]
 
 
                    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))
+#list-lookup should be faster than dict-lookup
 
 
-def PolygoniseTri(g, iso, v0, v1, v2, v3):
+def polygoniseTri(g, iso, v0, v1, v2, v3):
+    """
+    Polygonise a tetrahedron given its vertices within a cube
+    This is an alternative algorithm to polygonisegrid.
+    It results in a smoother surface but more triangular facets.
+
+                       + 0
+                      /|\
+                     / | \
+                    /  |  \
+                   /   |   \
+                  /    |    \
+                 /     |     \
+                +-------------+ 1
+               3 \     |     /
+                  \    |    /
+                   \   |   /
+                    \  |  /
+                     \ | /
+                      \|/
+                       + 2
+
+    It's main purpose is still to polygonise a gridded dataset and
+    would normally be called 6 times, one for each tetrahedron making
+    up the grid cell.
+    """
+
     #   Determine which of the 16 cases we have given which vertices
     #   are above or below the isosurface
-
     triindex = 0
-    if g.val[v0] < iso:
+    if g.val[v0] < iso:  # vertex 0 is "outside"
         triindex |= 1
-    if g.val[v1] < iso:
+    if g.val[v1] < iso:  # vertex 1 is "outside"
         triindex |= 2
-    if g.val[v2] < iso:
+    if g.val[v2] < iso:  # vertex 2 is "outside"
         triindex |= 4
-    if g.val[v3] < iso:
+    if g.val[v3] < iso:  # vertex 3 is "outside"
         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()
+    #TDOD: flip the triangle in some cases
     return l
 
 
 def make(density_function, isolevel, xmin, xmax, ymin, ymax, zmin, zmax):
+    """
+    Evaluate density_funktion on a regular grid and return the resulting triangles
+    triangles are generated for cells on the grid that are 
+    part denser than isolevel (inside) and part less dense then isolevel (outside)
+    """
     triangles = []
     p = [None] * 8
     val = [None] * 8
                 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))
+                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
 
+
+def make_simple(df, size, steps, isolevel):
+    stepscale = size / steps
+    dr = lambda x, y, z: df(x * stepscale, y * stepscale, z * stepscale)
+    sh = steps / 2
+    return make(dr, isolevel, -sh, sh, -sh, sh, -sh, sh)
+
+
+def make_from_cache(cached, isolevel):
+    #cached = cache(df, size, steps)
+    dr = lambda x, y, z: cached[x][y][z]
+    xmin = 0
+    xmax = len(cached) - 1
+    ymin = 0
+    ymax = len(cached[0]) - 1
+    zmin = 0
+    zmax = len(cached[0][0]) - 1
+    return make(dr, isolevel, xmin, xmax, ymin, ymax, zmin, zmax)
+
+
 if __name__ == "__main__":
-    data = readdata(Composite([Metaball(Vec3(0, 0, 0), 2), (lambda x, y, z:-Metaball(Vec3(1, 0, 0), 2)(x, y, z))]), 5, 41)
+    #this is a non-interactive example/test but you can (and should) just use standard controls
+
+    #some initial settings
+    df = Composite([Metaball(Vec3(0, 0, 0), 2), (lambda x, y, z: (-Metaball(Vec3(1, 0, 0), 2)(x, y, z)))])
     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)
+    size = 8.0
+    steps = 32
+    #generate triangles
+    tris = make_simple(df, size, steps, isolevel)
+
     import direct.directbase.DirectStart
-    tris = make(dr, isolevel, xmin, xmax, ymin, ymax, zmin, zmax)
+    #init meshdrawer as its the simplest way to draw from python
     md = MeshDrawer()
     mdNode = md.getRoot()
-    mdNode.setTwoSided(True)
+    mdNode.setTwoSided(True) #TODO: put the triangles the right direction on generation
     mdNode.reparentTo(render)  # @UndefinedVariable
     #mdNode.setRenderModeWireframe(True)
-    md.setBudget(xmax * ymax * zmax * 1) #actually 12 instead of 1 in general
+    md.setBudget(len(tris)) #we kow how many tris we have
     md.begin(base.cam, render)  # @UndefinedVariable
-    c = (1, 1, 1, 0.5)
-    uv = (0, 0)
+    c = (1, 1, 1, 0.5) #set a default color
+    uv = (0, 0) #we dont use textures in this example so tghis is there just because md.tri needs it
     for tri in tris:
         v1, v2, v3 = tri.p
         md.tri(v1, c, uv, v2, c, uv, v3, c, uv)