first pass at generalized dual contours

Author: sjkelly

REQUIRE

@@ -1,2 +1,6 @@
 julia 0.7
 GeometryTypes 0.4
+LinearAlgebra
+Roots
+StaticArrays
+ForwardDiff

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

@@ -1,22 +1,16 @@
-using LinearAlgebra
-using Roots
-using GeometryTypes
-using FileIO
-using StaticArrays
-using ForwardDiff
-
 #Cardinal directions
-dirs = [ [1,0,0], [0,1,0], [0,0,1] ]
+const dirs = ( Point(1,0,0), Point(0,1,0), Point(0,0,1) )
 
 #Vertices of cube
-cube_verts = [[0, 0, 0], [0, 0, 1], [0, 1, 0], [0, 1, 1],
-              [1, 0, 0], [1, 0, 1], [1, 1, 0], [1, 1, 1]]
+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
-cube_edges = [(0, 1), (0, 2), (0, 1), (0, 4), (0, 2), (0, 4), (2, 3), (1, 3),
+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)]
+              (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
@@ -32,24 +26,28 @@ end
 # f = implicit function
 # df = gradient of f
 # nc = resolution
-function dual_contour(f, nc)
+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, y in 0:nc, z in 0:nc
-        o = [x,y,z]
+    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)>0 for v in cube_verts ]
+        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], o+cube_verts[e[2]+1])
-            for e in cube_edges if cube_signs[e[1]+1] != cube_signs[e[2]+1] ]
+        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)
@@ -65,70 +63,44 @@ function dual_contour(f, nc)
         end
 
         #Emit one vertex per every cube that crosses
-        push!(vindex, o => length(dc_verts))
+        push!(vindex, idx => length(dc_verts))
         push!(dc_verts, (v, ForwardDiff.gradient(f,v)))
     end
 
     #Construct faces
-    dc_faces = []
-    for x in 0:nc, y in 0:nc, z in 0:nc
-        if !haskey(vindex,[x,y,z])
+    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 = [x,y,z]
+        o = Point(x,y,z) .* scale + orig
         for i in (1,2,3)
             for j in 1:i
-                if haskey(vindex,o + dirs[i]) && haskey(vindex,o + dirs[j]) && haskey(vindex,o + dirs[i] + dirs[j])
+                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[o]+1]
-                    v2, tn2 = dc_verts[vindex[o+dirs[i]]+1]
-                    v3, tn3 = dc_verts[vindex[o+dirs[j]]+1]
-                    @show v1,v2, v3
+                    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, [vindex[o], vindex[o+dirs[i]], vindex[o+dirs[j]]] )
-                        push!(dc_faces, [vindex[o+dirs[i]+dirs[j]], vindex[o+dirs[j]], vindex[o+dirs[i]]] )
+                        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, [vindex[o], vindex[o+dirs[j]], vindex[o+dirs[i]]] )
-                        push!(dc_faces, [vindex[o+dirs[i]+dirs[j]], vindex[o+dirs[i]], vindex[o+dirs[j]]] )
+                        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 dc_verts, dc_faces
-end
-
-
-center = [16,16,16]
-radius = 10
-
-function test_f(x)
-    d = x-center
-    return dot(d,d) - radius^2
-end
-
-function runion(x,y)
-    x + y - sqrt(x*x + y*y)
-end
-
-function softmax(x,y)
-    log(exp(3*x)+exp(3*y))/3
-end
-
-function test_sq(x)
-    d = x-center
-    softmax(softmax(-d[3], d[3]-5), d[1]*d[1] + d[2]*d[2] - radius )
+    return HomogenousMesh([Point(v[1]...) for v in dc_verts], dc_faces)
 end
 
-
-verts, tris = dual_contour(test_sq, 36)
-
-m = HomogenousMesh([Point(v[1]...) for v in verts], [Face(t[1]+1,t[2]+1,t[3]+1) for t in tris])
-
-save("test2.ply",m)

test/runtests.jl

@@ -105,6 +105,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))