Source

isosurfacesforp3d / src / tetrahedrons.py

Full commit
from panda3d.core import Vec3, MeshDrawer  # @UnresolvedImport
from density import Metaball, Composite
SMALL_DOUBLE = 0.00001


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):
    """
    A simple Triangle
    """
    def __init__(self, p1, p2, p3):
        self.p = [p1, p2, p3]  # vertices

    def flip(self):
        self.p.reverse()

    def __repr__(self):
        return "%s(%r, %r, %r)" % (self.__class__.__name__,
                                   self.p[0], self.p[1], self.p[2])


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 = []
    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 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) < SMALL_DOUBLE:
        return(p2)
    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]))
    ]


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))
#list-lookup should be faster than dict-lookup


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:  # vertex 0 is "outside"
        triindex |= 1
    if g.val[v1] < iso:  # vertex 1 is "outside"
        triindex |= 2
    if g.val[v2] < iso:  # vertex 2 is "outside"
        triindex |= 4
    if g.val[v3] < iso:  # vertex 3 is "outside"
        triindex |= 8

    l = trianglefs[triindex](g, iso, v0, v1, v2, v3)
    #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
    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


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__":
    #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
    size = 8.0
    steps = 32
    #generate triangles
    tris = make_simple(df, size, steps, isolevel)

    import direct.directbase.DirectStart
    #init meshdrawer as its the simplest way to draw from python
    md = MeshDrawer()
    mdNode = md.getRoot()
    mdNode.setTwoSided(True) #TODO: put the triangles the right direction on generation
    mdNode.reparentTo(render)  # @UndefinedVariable
    #mdNode.setRenderModeWireframe(True)
    md.setBudget(len(tris)) #we kow how many tris we have
    md.begin(base.cam, render)  # @UndefinedVariable
    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)
    md.end()
    run()  # @UndefinedVariable