From 4a389a30b00a74fe27febbf4d4f09bb683593eab Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Sat, 15 Jun 2019 18:38:18 -0400 Subject: [PATCH 1/5] wip dual contours --- Project.toml | 4 ++ src/Meshing.jl | 8 +++- src/dual_contours.jl | 106 +++++++++++++++++++++++++++++++++++++++++++ test/runtests.jl | 9 ++++ 4 files changed, 126 insertions(+), 1 deletion(-) create mode 100644 src/dual_contours.jl diff --git a/Project.toml b/Project.toml index afdc390..fd1851a 100644 --- a/Project.toml +++ b/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" diff --git a/src/Meshing.jl b/src/Meshing.jl index a8218a6..d0bfae6 100644 --- a/src/Meshing.jl +++ b/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 diff --git a/src/dual_contours.jl b/src/dual_contours.jl new file mode 100644 index 0000000..a8af7d0 --- /dev/null +++ b/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 + diff --git a/test/runtests.jl b/test/runtests.jl index e36573d..f2e2394 100644 --- a/test/runtests.jl +++ b/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)) From 027d167986971882e6636141a8c1bdac58fc5722 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Sun, 23 Jun 2019 16:19:33 -0400 Subject: [PATCH 2/5] add dual contours --- README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README.md b/README.md index 020e439..40f00d1 100644 --- a/README.md +++ b/README.md @@ -9,6 +9,7 @@ Algorithms included: * [Marching Tetrahedra](https://en.wikipedia.org/wiki/Marching_tetrahedra) * [Marching Cubes](https://en.wikipedia.org/wiki/Marching_cubes) * [Naive Surface Nets](https://0fps.net/2012/07/12/smooth-voxel-terrain-part-2/) +* [Dual Contours](https://en.wikipedia.org/wiki/Isosurface#Dual_Contouring) ## Interface @@ -44,6 +45,7 @@ Three meshing algorithms exist: * `MarchingCubes()` * `MarchingTetrahedra()` * `NaiveSurfaceNets()` +* `DualContours()` Each takes an optional `iso` and `eps` parameter, e.g. `MarchingCubes(0.0,1e-6)`. From 4fc63195f22f5b2eedcea06220a889521e05b06e Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Mon, 24 Jun 2019 20:56:27 -0400 Subject: [PATCH 3/5] integrate Dual contours into the API --- src/dual_contours.jl | 38 +++++++++++++++++++++++++++++--------- test/runtests.jl | 7 ++++--- 2 files changed, 33 insertions(+), 12 deletions(-) diff --git a/src/dual_contours.jl b/src/dual_contours.jl index a8af7d0..013c4e6 100644 --- a/src/dual_contours.jl +++ b/src/dual_contours.jl @@ -22,19 +22,21 @@ function estimate_hermite(f, v0, 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}) + +function dual_contours(f::Function, + bounds::HyperRectangle, + samples::NTuple{3,Int}=(128,128,128), + iso=0.0, + MT::Type{M}=SimpleMesh{Point{3,Float64},Face{3,Int}}, + eps=0.00001) where {ST,FT,M<:AbstractMesh} orig = origin(bounds) width = widths(bounds) - scale = width ./ Point(nc) + scale = width ./ Point(samples) #Compute vertices dc_verts = [] vindex = Dict() - for x in 0:nc[1], y in 0:nc[2], z in 0:nc[3] + for x in 0:samples[1], y in 0:samples[2], z in 0:samples[3] idx = Point(x,y,z) o = Point(x,y,z) .* scale + orig @@ -69,7 +71,7 @@ function dual_contours(f, bounds::HyperRectangle, nc::NTuple{3,Int}) #Construct faces dc_faces = Face[] - for x in 0:nc[1], y in 0:nc[2], z in 0:nc[3] + for x in 0:samples[1], y in 0:samples[2], z in 0:samples[3] idx = Point(x,y,z) if !haskey(vindex,idx) @@ -101,6 +103,24 @@ function dual_contours(f, bounds::HyperRectangle, nc::NTuple{3,Int}) end end - return HomogenousMesh([Point(v[1]...) for v in dc_verts], dc_faces) + return MT([Point(v[1]...) for v in dc_verts], dc_faces) +end + +struct DualContours{T} <: AbstractMeshingAlgorithm + iso::T + eps::T +end + +DualContours(iso::T1=0.0, eps::T2=1e-3) where {T1, T2} = DualContours{promote_type(T1, T2)}(iso, eps) + +function (::Type{MT})(df::SignedDistanceField, method::DualContours)::MT where {MT <: AbstractMesh} + dual_contours(df, method.iso, MT, method.eps) +end + +function (::Type{MT})(f::Function, h::HyperRectangle, size::NTuple{3,Int}, method::DualContours)::MT where {MT <: AbstractMesh} + dual_contours(f, h, size, method.iso, MT, method.eps) end +function (::Type{MT})(f::Function, h::HyperRectangle, method::DualContours; size::NTuple{3,Int}=(128,128,128))::MT where {MT <: AbstractMesh} + dual_contours(f, h, size, method.iso, MT, method.eps) +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index f2e2394..fbc3f67 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -145,9 +145,10 @@ using LinearAlgebra: dot, norm 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 + m = SimpleMesh(sphere, HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), (50,50,50), DualContours()) + # the current LLS approach in DualContours is not numerically robust and has varied outputs + @test_broken length(vertices(m)) == 11754 + @test_broken length(faces(m)) == 51186 end @testset "AbstractMeshingAlgorithm interface" begin From f0f21aed795973626787a58e7bcd8be8bb65b69f Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Tue, 25 Jun 2019 20:26:04 -0400 Subject: [PATCH 4/5] use pseudo inverse --- src/dual_contours.jl | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/dual_contours.jl b/src/dual_contours.jl index 013c4e6..125e1eb 100644 --- a/src/dual_contours.jl +++ b/src/dual_contours.jl @@ -57,7 +57,17 @@ function dual_contours(f::Function, A[i,:] = h_data[i][2] end b = [ dot(d[1],d[2]) for d in h_data ] - v = A\b + + # use SVD for pseudo inverse + # F = svd(A, full=true) + # S = zeros(size(A)) + # for i = 1:length(F.S) + # S[i,i] = 1/(F.S[i]) + # end + # ainv = F.V*S*F.U' + + #compute the vertex from SVD + v = pinv(A)*b #Throw out failed solutions if norm(v-o) > 2 From 90426adc62ec74eb48785ed2b30bfd63d504c945 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Tue, 25 Jun 2019 20:46:42 -0400 Subject: [PATCH 5/5] [test] fix dc test --- src/dual_contours.jl | 10 +--------- test/runtests.jl | 5 ++--- 2 files changed, 3 insertions(+), 12 deletions(-) diff --git a/src/dual_contours.jl b/src/dual_contours.jl index 125e1eb..88b42da 100644 --- a/src/dual_contours.jl +++ b/src/dual_contours.jl @@ -58,15 +58,7 @@ function dual_contours(f::Function, end b = [ dot(d[1],d[2]) for d in h_data ] - # use SVD for pseudo inverse - # F = svd(A, full=true) - # S = zeros(size(A)) - # for i = 1:length(F.S) - # S[i,i] = 1/(F.S[i]) - # end - # ainv = F.V*S*F.U' - - #compute the vertex from SVD + #compute the vertex using pseudo inverse v = pinv(A)*b #Throw out failed solutions diff --git a/test/runtests.jl b/test/runtests.jl index fbc3f67..e359e8a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -146,9 +146,8 @@ using LinearAlgebra: dot, norm sphere(v) = sqrt(sum(dot(v,v))) - 1 # sphere m = SimpleMesh(sphere, HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), (50,50,50), DualContours()) - # the current LLS approach in DualContours is not numerically robust and has varied outputs - @test_broken length(vertices(m)) == 11754 - @test_broken length(faces(m)) == 51186 + @test length(vertices(m)) == 11754 + @test length(faces(m)) == 51186 end @testset "AbstractMeshingAlgorithm interface" begin