working mesh sdf

examples/vdb.py

@@ -2,33 +2,51 @@
 import meshio
 import pyopenvdb as vdb
 import sys
+import threading
 
-def main():
-    args = sys.argv[1:]
-    if len(args) != 1:
-        print('Usage: python vdb.py input.stl')
-        return
-
-    mesh = meshio.read(args[0])
+@sdf3
+def mesh_from_file(path):
+    mesh = meshio.read(path)
 
     points = mesh.points
     triangles = mesh.cells[0].data
+    estimator = box(a=points.min(axis=0), b=points.max(axis=0))
+
+    lock = threading.Lock()
+    grids = {}
+
+    def get_grid(voxel_size):
+        with lock:
+            if voxel_size not in grids:
+                transform = vdb.createLinearTransform(voxelSize=voxel_size)
+                grid = vdb.FloatGrid.createLevelSetFromPolygons(
+                    points, triangles=triangles, transform=transform, halfWidth=3)
+                grids[voxel_size] = grid
+            return grids[voxel_size]
+
+    def f(p):
+        if not hasattr(p, 'info'):
+            return estimator(p)
+        nx = len(np.unique(p[:,0]))
+        ny = len(np.unique(p[:,1]))
+        nz = len(np.unique(p[:,2]))
+        voxel_size = float((p[1] - p[0]).max()) if len(p) > 1 else 1
+        voxel_size = round(voxel_size, 9)
+        transform = vdb.createLinearTransform(voxelSize=voxel_size)
+        grid = get_grid(voxel_size)
+        a = np.zeros(nx * ny * nz).reshape((nx, ny, nz))
+        ijk = transform.worldToIndex(p[0])
+        ijk = [int(round(x)) for x in ijk]
+        grid.copyToArray(a, ijk=ijk)
+        return a
+
+    return f
 
-    transform = vdb.createLinearTransform(voxelSize=0.25)
-
-    grid = vdb.FloatGrid.createLevelSetFromPolygons(points, triangles=triangles, transform=transform, halfWidth=3)
-
-    print(dir(grid))
-    print(grid.evalActiveVoxelBoundingBox())
-
-    points, triangles, quads = grid.convertToPolygons(adaptivity=0)
-    triangles = triangles[:,[2,1,0]]
-    triangles = np.concatenate([triangles, quads[:,[2,1,0]], quads[:,[0,3,2]]])
-    cells = [
-        ('triangle', triangles),
-    ]
-    mesh = meshio.Mesh(points, cells)
-    mesh.write('out.stl')
+def main():
+    f = mesh_from_file(sys.argv[1])
+    # f = f.shell(0.25)
+    f &= slab(y0=-2)
+    f.save('out.stl', step=0.2, sparse=False)
 
 if __name__ == '__main__':
     main()

sdf/mesh.py

@@ -13,6 +13,16 @@
 SAMPLES = 2 ** 22
 BATCH_SIZE = 32
 
+class InfoArray(np.ndarray):
+    def __new__(cls, input_array, info=None):
+        obj = np.asarray(input_array).view(cls)
+        obj.info = info
+        return obj
+    def __array_finalize__(self, obj):
+        if obj is None:
+            return
+        self.info = getattr(obj, 'info', None)
+
 def _marching_cubes(volume, level=0):
     verts, faces, _, _ = measure.marching_cubes(volume, level)
     return verts[faces].reshape((-1, 3))
@@ -48,6 +58,8 @@ def _worker(sdf, job, sparse):
         return None
         # return _debug_triangles(X, Y, Z)
     P = _cartesian_product(X, Y, Z)
+    # TODO: pass info (step, shape, bounds, etc.)
+    P = InfoArray(P, info=None)
     volume = sdf(P).reshape((len(X), len(Y), len(Z)))
     try:
         points = _marching_cubes(volume)