wip dual contours

Author: sjkelly

Project.toml

@@ -2,7 +2,11 @@ name = "Meshing"
 uuid = "e6723b4c-ebff-59f1-b4b7-d97aa5274f73"
 
 [deps]
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GeometryTypes = "4d00f742-c7ba-57c2-abde-4428a4b178cb"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [extras]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"

src/Meshing.jl

@@ -1,16 +1,22 @@
 module Meshing
 
 using GeometryTypes
+using Roots
+using ForwardDiff
+using LinearAlgebra
+using StaticArrays
 
 abstract type AbstractMeshingAlgorithm end
 
 include("marching_tetrahedra.jl")
 include("marching_cubes.jl")
 include("surface_nets.jl")
+include("dual_contours.jl")
 
 export marching_cubes,
        MarchingCubes,
        MarchingTetrahedra,
-       NaiveSurfaceNets
+       NaiveSurfaceNets,
+       DualContours
 
 end # module

src/dual_contours.jl

@@ -0,0 +1,106 @@
+#Cardinal directions
+const dirs = ( Point(1,0,0), Point(0,1,0), Point(0,0,1) )
+
+#Vertices of cube
+const cube_verts = (Point(0, 0, 0), Point(0, 0, 1), Point(0, 1, 0),
+                    Point(0, 1, 1), Point(1, 0, 0), Point(1, 0, 1),
+                    Point(1, 1, 0), Point(1, 1, 1))
+
+
+#Edges of cube
+const cube_edges_dc = ((0, 1), (0, 2), (0, 1), (0, 4), (0, 2), (0, 4), (2, 3), (1, 3),
+              (4, 5), (1, 5), (4, 6), (2, 6), (4, 5), (4, 6), (2, 3), (2, 6),
+              (1, 3), (1, 5), (6, 7), (5, 7), (6, 7), (3, 7), (5, 7), (3, 7))
+
+
+#Use non-linear root finding to compute intersection point
+function estimate_hermite(f, v0, v1)
+    l(t) = f((1.0-t)*v0 + t*v1)
+    dl(t) = ForwardDiff.derivative(l,t)
+    t0 = find_zero((l,dl),(0, 1))
+    x0 = (1.0-t0)*v0 + t0*v1
+    return (x0, ForwardDiff.gradient(f,x0))
+end
+
+#Input:
+# f = implicit function
+# df = gradient of f
+# nc = resolution
+function dual_contours(f, bounds::HyperRectangle, nc::NTuple{3,Int})
+
+    orig = origin(bounds)
+    width = widths(bounds)
+    scale = width ./ Point(nc)
+    #Compute vertices
+    dc_verts = []
+    vindex   = Dict()
+    for x in 0:nc[1], y in 0:nc[2], z in 0:nc[3]
+        idx = Point(x,y,z)
+        o = Point(x,y,z) .* scale + orig
+
+        #Get signs for cube
+        cube_signs = [ f(o+(v.*scale))>0 for v in cube_verts ]
+
+        if all(cube_signs) || !any(cube_signs)
+            continue
+        end
+
+        #Estimate hermite data
+        h_data = [ estimate_hermite(f, o+cube_verts[e[1]+1].*scale, o+cube_verts[e[2]+1].*scale)
+            for e in cube_edges_dc if cube_signs[e[1]+1] != cube_signs[e[2]+1] ]
+
+        #Solve qef to get vertex
+        A = Array{Float64,2}(undef,length(h_data),3)
+        for i in eachindex(h_data)
+            A[i,:] = h_data[i][2]
+        end
+        b = [ dot(d[1],d[2]) for d in h_data ]
+        v = A\b
+
+        #Throw out failed solutions
+        if norm(v-o) > 2
+            continue
+        end
+
+        #Emit one vertex per every cube that crosses
+        push!(vindex, idx => length(dc_verts))
+        push!(dc_verts, (v, ForwardDiff.gradient(f,v)))
+    end
+
+    #Construct faces
+    dc_faces = Face[]
+    for x in 0:nc[1], y in 0:nc[2], z in 0:nc[3]
+
+        idx = Point(x,y,z)
+        if !haskey(vindex,idx)
+            continue
+        end
+
+        #Emit one face per each edge that crosses
+        o = Point(x,y,z) .* scale + orig
+        for i in (1,2,3)
+            for j in 1:i
+                if haskey(vindex,idx+dirs[i]) && haskey(vindex,idx + dirs[j]) && haskey(vindex,idx + dirs[i] + dirs[j])
+                    # determine orientation of the face from the true normal
+                    v1, tn1 = dc_verts[vindex[idx]+1]
+                    v2, tn2 = dc_verts[vindex[idx+dirs[i]]+1]
+                    v3, tn3 = dc_verts[vindex[idx+dirs[j]]+1]
+
+                    e1 = v1-v2
+                    e2 = v1-v3
+                    c = cross(e1,e2)
+                    if dot(c, tn1) > 0
+                        push!(dc_faces, Face{3,Int}(vindex[idx]+1, vindex[idx+dirs[i]]+1, vindex[idx+dirs[j]]+1) )
+                        push!(dc_faces, Face{3,Int}(vindex[idx+dirs[i]+dirs[j]]+1, vindex[idx+dirs[j]]+1, vindex[idx+dirs[i]]+1) )
+                    else
+                        push!(dc_faces, Face{3,Int}(vindex[idx]+1, vindex[idx+dirs[j]]+1, vindex[idx+dirs[i]]+1) )
+                        push!(dc_faces, Face{3,Int}(vindex[idx+dirs[i]+dirs[j]]+1, vindex[idx+dirs[i]]+1, vindex[idx+dirs[j]]+1) )
+                    end
+                end
+            end
+        end
+
+    end
+    return HomogenousMesh([Point(v[1]...) for v in dc_verts], dc_faces)
+end
+

test/runtests.jl

@@ -141,6 +141,15 @@ using LinearAlgebra: dot, norm
         @test minimum(vertices(mesh)) ≈ [-0.5, -0.5, -0.5] atol=0.2
     end
 
+    @testset "Dual Contours" begin
+
+        sphere(v) = sqrt(sum(dot(v,v))) - 1 # sphere
+
+        m = Meshing.dual_contours(sphere, HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), (50,50,50))
+        @test length(vertices(m)) == 11754
+        @test length(faces(m)) == 51186
+    end
+
     @testset "AbstractMeshingAlgorithm interface" begin
         f = x -> norm(x) - 0.5
         bounds = HyperRectangle(Vec(-1, -1, -1), Vec(2, 2, 2))