[WIP] Adaptive Distance Field meshing

Project.toml

@@ -2,7 +2,10 @@ name = "Meshing"
 uuid = "e6723b4c-ebff-59f1-b4b7-d97aa5274f73"
 
 [deps]
+AdaptiveDistanceFields = "a1957575-6125-5dba-8f92-417d2d1f4a46"
 GeometryTypes = "4d00f742-c7ba-57c2-abde-4428a4b178cb"
+RegionTrees = "dee08c22-ab7f-5625-9660-a9af2021b33f"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [extras]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"

src/Meshing.jl

@@ -1,14 +1,18 @@
 module Meshing
 
 using GeometryTypes
+using AdaptiveDistanceFields
+import RegionTrees
+using StaticArrays
 
 abstract type AbstractMeshingAlgorithm end
 
 include("marching_tetrahedra.jl")
 include("marching_cubes.jl")
 include("surface_nets.jl")
 
-export marching_cubes,
+export AdaptiveMarchingCubes,
+       marching_cubes,
        MarchingCubes,
        MarchingTetrahedra,
        NaiveSurfaceNets

src/marching_cubes.jl

@@ -443,6 +443,71 @@ function marching_cubes(f::Function,
     MT(vts,fcs)
 end
 
+function marching_cubes_adf(f::Function,
+                        origin,
+                        widths,
+                        rtol=1e-2,
+                        atol=1e-2,
+                        iso=0.0,
+                        MT::Type{M}=SimpleMesh{Point{3,Float64},Face{3,Int}},
+                        eps=0.00001) where {ST,FT,M<:AbstractMesh}
+
+    adf = AdaptiveDistanceField(f, origin, widths, rtol, atol)
+
+    # arrays for vertices and faces
+    vts = Point{3,Float64}[]
+    fcs = Face{3,Int}[]
+    vertlist = Vector{Point{3,Float64}}(undef, 12)
+    leaves = RegionTrees.allleaves(adf.root)
+
+    @inbounds for leaf in leaves
+        o = leaf.boundary.origin
+        w = leaf.boundary.widths
+        points = (Point{3,Float64}(o),
+                  o + w .* Point{3,Float64}(1,0,0),
+                  o + w .* Point{3,Float64}(1,1,0),
+                  o + w .* Point{3,Float64}(0,1,0),
+                  o + w .* Point{3,Float64}(0,0,1),
+                  o + w .* Point{3,Float64}(1,0,1),
+                  o + w .* Point{3,Float64}(1,1,1),
+                  o + w .* Point{3,Float64}(0,1,1))
+
+        iso_vals = map(adf,points)
+
+        #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
+        find_vertices_interp!(vertlist, points, iso_vals, cubeindex, iso, eps)
+
+        # 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
+
 @inline function find_vertices_interp!(vertlist, points, iso_vals, cubeindex, iso, eps)
      if (edge_table[cubeindex] & 1 != 0)
      vertlist[1] =
@@ -514,6 +579,19 @@ end
 
 MarchingCubes(iso::T1=0.0, eps::T2=1e-3) where {T1, T2} = MarchingCubes{promote_type(T1, T2)}(iso, eps)
 
+struct AdaptiveMarchingCubes{T} <: AbstractMeshingAlgorithm
+     iso::T
+     eps::T
+     rtol::T
+     atol::T
+end
+
+AdaptiveMarchingCubes(iso::T1=0.0, eps::T2=1e-3, rtol::T3=1e-2, atol::T4=1e-2) where {T1, T2, T3, T4} = AdaptiveMarchingCubes{promote_type(T1, T2, T3, T4)}(iso, eps, rtol, atol)
+
+#
+# Marching Cubes
+#
+
 function (::Type{MT})(df::SignedDistanceField, method::MarchingCubes)::MT where {MT <: AbstractMesh}
      marching_cubes(df, method.iso, MT, method.eps)
 end
@@ -524,4 +602,16 @@ end
 
 function (::Type{MT})(f::Function, h::HyperRectangle, method::MarchingCubes; size::NTuple{3,Int}=(128,128,128))::MT where {MT <: AbstractMesh}
      marching_cubes(f, h, size, method.iso, MT, method.eps)
+end
+
+#
+# Adaptive Marching Cubes
+#
+
+function (::Type{MT})(f::Function, h::HyperRectangle, method::AdaptiveMarchingCubes)::MT where {MT <: AbstractMesh}
+     marching_cubes_adf(f, SVector(h.origin), SVector(h.widths), method.rtol, method.atol, method.iso, MT, method.eps)
+end
+
+function (::Type{MT})(f::Function, h::HyperRectangle, samples::NTuple, method::AdaptiveMarchingCubes)::MT where {MT <: AbstractMesh}
+     marching_cubes_adf(f, SVector(h.origin), SVector(h.widths), method.rtol, method.atol, method.iso, MT, method.eps)
 end

test/runtests.jl

@@ -112,6 +112,16 @@ using LinearAlgebra: dot, norm
         @test length(faces(m)) == length(faces(mf))
     end
 
+    @testset "adaptive marching cubes" begin
+        m = SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), AdaptiveMarchingCubes()) do v
+            sqrt(sum(dot(v,v))) - 1 # sphere
+        end
+
+        @show length(vertices(m))
+        @show length(faces(m))
+
+    end
+
     @testset "respect origin" begin
         # verify that when we construct a mesh, that mesh:
         #   a) respects the origin of the SDF