[MC] First pass at allocation-free meshing for MC

Author: sjkelly

src/marching_cubes.jl

@@ -392,6 +392,144 @@ function marching_cubes(sdf::SignedDistanceField{3,ST,FT},
     MT(vts,fcs)
 end
 
+
+function marching_cubes(f::Function,
+                        bounds::HyperRectangle,
+                        samples::NTuple{3,Int}=(256,256,256),
+                        iso=0.0,
+                        MT::Type{M}=SimpleMesh{Point{3,Float64},Face{3,Int}},
+                        eps=0.00001) where {ST,FT,M<:AbstractMesh}
+    nx, ny, nz = samples[1], samples[2], samples[3]
+    w = widths(bounds)
+    orig = origin(bounds)
+
+    # we subtract one from the length along each axis because
+    # an NxNxN SDF has N-1 cells on each axis
+    s = Point{3,Float64}(w[1]/(nx-1), w[2]/(ny-1), w[3]/(nz-1))
+
+    # arrays for vertices and faces
+    vts = Point{3,Float64}[]
+    fcs = Face{3,Int}[]
+    mt = max(nx,ny,nz)
+    sizehint!(vts, mt*mt*6)
+    sizehint!(fcs, mt*mt*2)
+    vertlist = Vector{Point{3,Float64}}(undef, 12)
+    iso_vals = Vector{Float64}(undef,8)
+    points = Vector{Point{3,Float64}}(undef,8)
+    @inbounds for xi = 1:nx-1, yi = 1:ny-1, zi = 1:nz-1
+
+
+        if zi == 1
+            points[1] = Point{3,Float64}(xi-1,yi-1,zi-1) .* s .+ orig
+            points[2] = Point{3,Float64}(xi,yi-1,zi-1) .* s .+ orig
+            points[3] = Point{3,Float64}(xi,yi,zi-1) .* s .+ orig
+            points[4] = Point{3,Float64}(xi-1,yi,zi-1) .* s .+ orig
+            points[5] = Point{3,Float64}(xi-1,yi-1,zi) .* s .+ orig
+            points[6] = Point{3,Float64}(xi,yi-1,zi) .* s .+ orig
+            points[7] = Point{3,Float64}(xi,yi,zi) .* s .+ orig
+            points[8] = Point{3,Float64}(xi-1,yi,zi) .* s .+ orig
+            for i = 1:8
+                iso_vals[i] = f(points[i])
+            end
+        else
+            points[1] = points[5]
+            points[2] = points[6]
+            points[3] = points[7]
+            points[4] = points[8]
+            points[5] = Point{3,Float64}(xi-1,yi-1,zi) .* s .+ orig
+            points[6] = Point{3,Float64}(xi,yi-1,zi) .* s .+ orig
+            points[7] = Point{3,Float64}(xi,yi,zi) .* s .+ orig
+            points[8] = Point{3,Float64}(xi-1,yi,zi) .* s .+ orig
+            iso_vals[1] = iso_vals[5]
+            iso_vals[2] = iso_vals[6]
+            iso_vals[3] = iso_vals[7]
+            iso_vals[4] = iso_vals[8]
+            iso_vals[5] = f(points[5])
+            iso_vals[6] = f(points[6])
+            iso_vals[7] = f(points[7])
+            iso_vals[8] = f(points[8])
+        end
+
+        #Determine the index into the edge table which
+        #tells us which vertices are inside of the surface
+        cubeindex = iso_vals[1] < iso ? 1 : 0
+        iso_vals[2] < iso && (cubeindex |= 2)
+        iso_vals[3] < iso && (cubeindex |= 4)
+        iso_vals[4] < iso && (cubeindex |= 8)
+        iso_vals[5] < iso && (cubeindex |= 16)
+        iso_vals[6] < iso && (cubeindex |= 32)
+        iso_vals[7] < iso && (cubeindex |= 64)
+        iso_vals[8] < iso && (cubeindex |= 128)
+        cubeindex += 1
+
+        # Cube is entirely in/out of the surface
+        edge_table[cubeindex] == 0 && continue
+
+        # Find the vertices where the surface intersects the cube
+        # TODO this can use the underlying function to find the zero.
+        # The underlying space is non-linear so there will be error otherwise
+        if (edge_table[cubeindex] & 1 != 0)
+          vertlist[1] =
+             vertex_interp(iso,points[1],points[2],iso_vals[1],iso_vals[2], eps)
+        end
+        if (edge_table[cubeindex] & 2 != 0)
+          vertlist[2] =
+             vertex_interp(iso,points[2],points[3],iso_vals[2],iso_vals[3], eps)
+        end
+        if (edge_table[cubeindex] & 4 != 0)
+          vertlist[3] =
+             vertex_interp(iso,points[3],points[4],iso_vals[3],iso_vals[4], eps)
+        end
+        if (edge_table[cubeindex] & 8 != 0)
+          vertlist[4] =
+             vertex_interp(iso,points[4],points[1],iso_vals[4],iso_vals[1], eps)
+        end
+        if (edge_table[cubeindex] & 16 != 0)
+          vertlist[5] =
+             vertex_interp(iso,points[5],points[6],iso_vals[5],iso_vals[6], eps)
+        end
+        if (edge_table[cubeindex] & 32 != 0)
+          vertlist[6] =
+             vertex_interp(iso,points[6],points[7],iso_vals[6],iso_vals[7], eps)
+        end
+        if (edge_table[cubeindex] & 64 != 0)
+          vertlist[7] =
+             vertex_interp(iso,points[7],points[8],iso_vals[7],iso_vals[8], eps)
+        end
+        if (edge_table[cubeindex] & 128 != 0)
+          vertlist[8] =
+             vertex_interp(iso,points[8],points[5],iso_vals[8],iso_vals[5], eps)
+        end
+        if (edge_table[cubeindex] & 256 != 0)
+          vertlist[9] =
+             vertex_interp(iso,points[1],points[5],iso_vals[1],iso_vals[5], eps)
+        end
+        if (edge_table[cubeindex] & 512 != 0)
+          vertlist[10] =
+             vertex_interp(iso,points[2],points[6],iso_vals[2],iso_vals[6], eps)
+        end
+        if (edge_table[cubeindex] & 1024 != 0)
+          vertlist[11] =
+             vertex_interp(iso,points[3],points[7],iso_vals[3],iso_vals[7], eps)
+        end
+        if (edge_table[cubeindex] & 2048 != 0)
+          vertlist[12] =
+             vertex_interp(iso,points[4],points[8],iso_vals[4],iso_vals[8], eps)
+        end
+
+        # Create the triangle
+        for i = 1:3:13
+            tri_table[cubeindex][i] == -1 && break
+            push!(vts, vertlist[tri_table[cubeindex][i  ]])
+            push!(vts, vertlist[tri_table[cubeindex][i+1]])
+            push!(vts, vertlist[tri_table[cubeindex][i+2]])
+            fct = length(vts)
+            push!(fcs, Face{3,Int}(fct, fct-1, fct-2))
+        end
+    end
+    MT(vts,fcs)
+end
+
 # Linearly interpolate the position where an isosurface cuts
 # an edge between two vertices, each with their own scalar value
 function vertex_interp(iso, p1, p2, valp1, valp2, eps = 0.00001)
@@ -415,3 +553,7 @@ MarchingCubes(iso::T1=0.0, eps::T2=1e-3) where {T1, T2} = MarchingCubes{promote_
 function (::Type{MT})(df::SignedDistanceField, method::MarchingCubes)::MT where {MT <: AbstractMesh}
      marching_cubes(df, method.iso, MT, method.eps)
 end
+
+function (::Type{MT})(f::Function, h::HyperRectangle, size::NTuple{3,Int}, method::MarchingCubes)::MT where {MT <: AbstractMesh}
+     marching_cubes(f, h, size, method.iso, MT, method.eps)
+end

test/runtests.jl

@@ -76,6 +76,15 @@ using LinearAlgebra: dot, norm
         @test m == m2
 
     end
+
+    @testset "marching cubes function" begin
+        m = marching_cubes(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)),(21,21,21)) do v
+            sqrt(sum(dot(v,v))) - 1 # sphere
+        end
+        @test length(vertices(m)) == 10968
+        @test length(faces(m)) == 3656
+    end
+
     @testset "respect origin" begin
         # verify that when we construct a mesh, that mesh:
         #   a) respects the origin of the SDF