# Commits

committed 1190039

• Participants
• Parent commits c1c4d2b

# 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)`