Source

isosurfacesforp3d / src / tetrahedrons.py

Full commit
from panda3d.core import Vec3, MeshDrawer # @UnresolvedImport

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 __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()

    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):
        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):
    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):
    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):
    #   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 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__":
    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