diff --git a/.travis.yml b/.travis.yml index def1521..01079fb 100644 --- a/.travis.yml +++ b/.travis.yml @@ -18,4 +18,4 @@ script: - if [[ -a .git/shallow ]]; then git fetch --unshallow; fi - julia -e 'Pkg.clone(pwd()); Pkg.build("Meshing"); Pkg.test("Meshing"; coverage=true)' after_success: - - julia -e 'cd(Pkg.dir("Meshing")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(process_folder()); Codecov.submit(process_folder())' + - julia -e 'cd(Pkg.dir("Meshing")); Pkg.add("Coverage"); using Coverage; Codecov.submit(process_folder())' diff --git a/src/Meshing.jl b/src/Meshing.jl index 0e158a0..017aa50 100644 --- a/src/Meshing.jl +++ b/src/Meshing.jl @@ -4,9 +4,13 @@ module Meshing using GeometryTypes +abstract type AbstractMeshingAlgorithm end + include("marching_tetrahedra.jl") include("marching_cubes.jl") -export marching_cubes +export marching_cubes, + MarchingCubes, + MarchingTetrahedra end # module diff --git a/src/marching_cubes.jl b/src/marching_cubes.jl index 1df0d9c..a6da327 100644 --- a/src/marching_cubes.jl +++ b/src/marching_cubes.jl @@ -287,12 +287,17 @@ However it may generate non-manifold meshes, while Marching Tetrahedra guarentees a manifold mesh. """ function marching_cubes(sdf::SignedDistanceField{3,ST,FT}, - iso=0.0, - MT::Type{M}=SimpleMesh{Point{3,Float64},Face{3,Int}}) where {ST,FT,M<:AbstractMesh} + iso=0.0, + MT::Type{M}=SimpleMesh{Point{3,Float64},Face{3,Int}}, + eps=0.00001) where {ST,FT,M<:AbstractMesh} nx, ny, nz = size(sdf) h = HyperRectangle(sdf) w = widths(h) - s = Point{3,Float64}(w[1]/nx, w[2]/ny, w[3]/nz) + orig = origin(HyperRectangle(sdf)) + + # 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}[] @@ -315,63 +320,63 @@ function marching_cubes(sdf::SignedDistanceField{3,ST,FT}, # Cube is entirely in/out of the surface edge_table[cubeindex] == 0 && continue - points = (Point{3,Float64}(xi-1,yi-1,zi-1).*s, - Point{3,Float64}(xi,yi-1,zi-1).*s, - Point{3,Float64}(xi,yi,zi-1).*s, - Point{3,Float64}(xi-1,yi,zi-1).*s, - Point{3,Float64}(xi-1,yi-1,zi).*s, - Point{3,Float64}(xi,yi-1,zi).*s, - Point{3,Float64}(xi,yi,zi).*s, - Point{3,Float64}(xi-1,yi,zi).*s) + points = (Point{3,Float64}(xi-1,yi-1,zi-1) .* s .+ orig, + Point{3,Float64}(xi,yi-1,zi-1) .* s .+ orig, + Point{3,Float64}(xi,yi,zi-1) .* s .+ orig, + Point{3,Float64}(xi-1,yi,zi-1) .* s .+ orig, + Point{3,Float64}(xi-1,yi-1,zi) .* s .+ orig, + Point{3,Float64}(xi,yi-1,zi) .* s .+ orig, + Point{3,Float64}(xi,yi,zi) .* s .+ orig, + Point{3,Float64}(xi-1,yi,zi) .* s .+ orig) # Find the vertices where the surface intersects the cube if (edge_table[cubeindex] & 1 != 0) vertlist[1] = - vertex_interp(iso,points[1],points[2],sdf[xi,yi,zi],sdf[xi+1,yi,zi]) + vertex_interp(iso,points[1],points[2],sdf[xi,yi,zi],sdf[xi+1,yi,zi], eps) end if (edge_table[cubeindex] & 2 != 0) vertlist[2] = - vertex_interp(iso,points[2],points[3],sdf[xi+1,yi,zi],sdf[xi+1,yi+1,zi]) + vertex_interp(iso,points[2],points[3],sdf[xi+1,yi,zi],sdf[xi+1,yi+1,zi], eps) end if (edge_table[cubeindex] & 4 != 0) vertlist[3] = - vertex_interp(iso,points[3],points[4],sdf[xi+1,yi+1,zi],sdf[xi,yi+1,zi]) + vertex_interp(iso,points[3],points[4],sdf[xi+1,yi+1,zi],sdf[xi,yi+1,zi], eps) end if (edge_table[cubeindex] & 8 != 0) vertlist[4] = - vertex_interp(iso,points[4],points[1],sdf[xi,yi+1,zi],sdf[xi,yi,zi]) + vertex_interp(iso,points[4],points[1],sdf[xi,yi+1,zi],sdf[xi,yi,zi], eps) end if (edge_table[cubeindex] & 16 != 0) vertlist[5] = - vertex_interp(iso,points[5],points[6],sdf[xi,yi,zi+1],sdf[xi+1,yi,zi+1]) + vertex_interp(iso,points[5],points[6],sdf[xi,yi,zi+1],sdf[xi+1,yi,zi+1], eps) end if (edge_table[cubeindex] & 32 != 0) vertlist[6] = - vertex_interp(iso,points[6],points[7],sdf[xi+1,yi,zi+1],sdf[xi+1,yi+1,zi+1]) + vertex_interp(iso,points[6],points[7],sdf[xi+1,yi,zi+1],sdf[xi+1,yi+1,zi+1], eps) end if (edge_table[cubeindex] & 64 != 0) vertlist[7] = - vertex_interp(iso,points[7],points[8],sdf[xi+1,yi+1,zi+1],sdf[xi,yi+1,zi+1]) + vertex_interp(iso,points[7],points[8],sdf[xi+1,yi+1,zi+1],sdf[xi,yi+1,zi+1], eps) end if (edge_table[cubeindex] & 128 != 0) vertlist[8] = - vertex_interp(iso,points[8],points[5],sdf[xi,yi+1,zi+1],sdf[xi,yi,zi+1]) + vertex_interp(iso,points[8],points[5],sdf[xi,yi+1,zi+1],sdf[xi,yi,zi+1], eps) end if (edge_table[cubeindex] & 256 != 0) vertlist[9] = - vertex_interp(iso,points[1],points[5],sdf[xi,yi,zi],sdf[xi,yi,zi+1]) + vertex_interp(iso,points[1],points[5],sdf[xi,yi,zi],sdf[xi,yi,zi+1], eps) end if (edge_table[cubeindex] & 512 != 0) vertlist[10] = - vertex_interp(iso,points[2],points[6],sdf[xi+1,yi,zi],sdf[xi+1,yi,zi+1]) + vertex_interp(iso,points[2],points[6],sdf[xi+1,yi,zi],sdf[xi+1,yi,zi+1], eps) end if (edge_table[cubeindex] & 1024 != 0) vertlist[11] = - vertex_interp(iso,points[3],points[7],sdf[xi+1,yi+1,zi],sdf[xi+1,yi+1,zi+1]) + vertex_interp(iso,points[3],points[7],sdf[xi+1,yi+1,zi],sdf[xi+1,yi+1,zi+1], eps) end if (edge_table[cubeindex] & 2048 != 0) vertlist[12] = - vertex_interp(iso,points[4],points[8],sdf[xi,yi+1,zi],sdf[xi,yi+1,zi+1]) + vertex_interp(iso,points[4],points[8],sdf[xi,yi+1,zi],sdf[xi,yi+1,zi+1], eps) end # Create the triangle @@ -399,3 +404,14 @@ function vertex_interp(iso, p1, p2, valp1, valp2, eps = 0.00001) return p end + +struct MarchingCubes{T} <: AbstractMeshingAlgorithm + iso::T + eps::T +end + +MarchingCubes(iso::T1=0.0, eps::T2=1e-3) where {T1, T2} = MarchingCubes{promote_type(T1, T2)}(iso, eps) + +function (::Type{MT})(df::SignedDistanceField, method::MarchingCubes)::MT where {MT <: AbstractMesh} + marching_cubes(df, method.iso, MT, method.eps) +end diff --git a/src/marching_tetrahedra.jl b/src/marching_tetrahedra.jl index ae9cc0c..e831f92 100644 --- a/src/marching_tetrahedra.jl +++ b/src/marching_tetrahedra.jl @@ -292,8 +292,24 @@ end isosurface(lsf,isoval) = isosurface(lsf,isoval, convert(eltype(lsf), 0.001)) +""" +The marchingTetrahedra function returns vertices on the (1-based) indices of the +SDF's data, ignoring its actual bounds. This function adjusts the vertices in +place so that they correspond to points within the SDF bounds. +""" +function _correct_vertices!(vts, sdf::SignedDistanceField) + bounds = HyperRectangle(sdf) + orig = origin(bounds) + w = widths(bounds) + s = w ./ Point(size(sdf) .- 1) # subtract 1 because an SDF with N points per side has N-1 cells + for i in eachindex(vts) + vts[i] = (vts[i] .- 1) .* s .+ orig # subtract 1 to fix 1-indexing + end +end + function (::Type{MT})(volume::Array{T, 3}, iso_val::Real, eps_val=0.001) where {MT <: AbstractMesh, T} + Base.depwarn("(mesh type)(volume::Array, iso_val, eps_val) is deprecated. Please use: (mesh type)(volume, MarchingTetrahedra(iso_val, eps_val))", :Meshing_MT_array_eps) iso_val = convert(T, iso_val) eps_val = convert(T, eps_val) vts, fcs = isosurface(volume, iso_val, eps_val) @@ -301,7 +317,26 @@ function (::Type{MT})(volume::Array{T, 3}, iso_val::Real, eps_val=0.001) where { end function (::Type{MT})(df::SignedDistanceField, eps_val=0.001) where MT <: AbstractMesh + Base.depwarn("(mesh type)(sdf::SignedDistanceField, eps_val) is deprecated. Please use: (mesh type)(sdf, MarchingTetrahedra(0, eps))", :Meshing_MT_sdf_eps) vts, fcs = isosurface(df.data, 0.0, eps_val) + _correct_vertices!(vts, df) MT(vts, fcs) end +struct MarchingTetrahedra{T} <: AbstractMeshingAlgorithm + iso::T + eps::T +end + +MarchingTetrahedra(iso::T1=0.0, eps::T2=1e-3) where {T1, T2} = MarchingTetrahedra{promote_type(T1, T2)}(iso, eps) + +function (::Type{MT})(sdf::SignedDistanceField, method::MarchingTetrahedra)::MT where {MT <: AbstractMesh} + vts, fcs = isosurface(sdf.data, method.iso, method.eps) + _correct_vertices!(vts, sdf) + MT(vts, fcs) +end + +function (::Type{MT}){MT <: AbstractMesh, T}(volume::Array{T, 3}, method::MarchingTetrahedra)::MT + vts, fcs = isosurface(volume, convert(T, method.iso), convert(T, method.eps)) + MT(vts, fcs) +end diff --git a/test/runtests.jl b/test/runtests.jl index d587201..5d34e5c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,46 +2,100 @@ using Meshing using Base.Test using GeometryTypes -# Produce a level set function that is a noisy version of the distance from -# the origin (such that level sets are noisy spheres). -# -# The noise should exercise marching tetrahedra's ability to produce a water- -# tight surface in all cases (unlike standard marching cubes). -# -N = 10 -sigma = 1.0 -distance = Float32[ sqrt(Float32(i*i+j*j+k*k)) for i = -N:N, j = -N:N, k = -N:N ] -distance = distance + sigma*rand(2*N+1,2*N+1,2*N+1) - -# Extract an isosurface. -# -lambda = N-2*sigma # isovalue - -msh = HomogenousMesh(distance,lambda) - -s2 = SignedDistanceField(HyperRectangle(Vec(0,0,0.),Vec(1,1,1.))) do v - sqrt(sum(dot(v,v))) - 1 # sphere -end -msh = HomogenousMesh(s2) -@test length(vertices(msh)) == 973 -@test length(faces(msh)) == 1830 +@testset "meshing" begin + @testset "noisy spheres" begin + # Produce a level set function that is a noisy version of the distance from + # the origin (such that level sets are noisy spheres). + # + # The noise should exercise marching tetrahedra's ability to produce a water- + # tight surface in all cases (unlike standard marching cubes). + # + N = 10 + sigma = 1.0 + distance = Float32[ sqrt(Float32(i*i+j*j+k*k)) for i = -N:N, j = -N:N, k = -N:N ] + distance = distance + sigma*rand(2*N+1,2*N+1,2*N+1) + + # Extract an isosurface. + # + lambda = N-2*sigma # isovalue + + msh = HomogenousMesh(distance, MarchingTetrahedra(lambda)) + + s2 = SignedDistanceField(HyperRectangle(Vec(0,0,0.),Vec(1,1,1.))) do v + sqrt(sum(dot(v,v))) - 1 # sphere + end + + msh = HomogenousMesh(s2, MarchingTetrahedra()) + @test length(vertices(msh)) == 973 + @test length(faces(msh)) == 1830 + end + + @testset "vertex interpolation" begin + @test Meshing.vertex_interp(0, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,0.5,0) + @test Meshing.vertex_interp(-1, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,0,0) + @test Meshing.vertex_interp(1, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,1,0) + end + + @testset "marching cubes" begin + sdf = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v + sqrt(sum(dot(v,v))) - 1 # sphere + end + + m = marching_cubes(sdf,0) + m2 = marching_cubes(sdf) + @test length(vertices(m)) == 10968 + @test length(faces(m)) == 3656 + @test m == m2 + + end + @testset "respect origin" begin + # verify that when we construct a mesh, that mesh: + # a) respects the origin of the SDF + # b) is correctly spaced so that symmetric functions create symmetric meshes + f = x -> norm(x) + bounds = HyperRectangle(Vec(-1, -1, -1), Vec(2, 2, 2)) + resolution = 0.1 + sdf = SignedDistanceField(f, bounds, resolution) + + for algorithm in (MarchingCubes(0.5), MarchingTetrahedra(0.5)) + mesh = @inferred GLNormalMesh(sdf, algorithm) + # should be centered on the origin + @test mean(vertices(mesh)) ≈ [0, 0, 0] atol=0.15*resolution + # and should be symmetric about the origin + @test maximum(vertices(mesh)) ≈ [0.5, 0.5, 0.5] + @test minimum(vertices(mesh)) ≈ [-0.5, -0.5, -0.5] + end + end + + @testset "AbstractMeshingAlgorithm interface" begin + f = x -> norm(x) - 0.5 + bounds = HyperRectangle(Vec(-1, -1, -1), Vec(2, 2, 2)) + resolution = 0.1 + sdf = SignedDistanceField(f, bounds, resolution) + + @testset "marching cubes" begin + m1 = @test_nowarn marching_cubes(sdf, 0.0, GLNormalMesh) + m2 = @test_nowarn GLNormalMesh(sdf, MarchingCubes()) + @test vertices(m1) == vertices(m2) + @test faces(m1) == faces(m2) + end -# Vertex interpolation -@test Meshing.vertex_interp(0, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,0.5,0) -@test Meshing.vertex_interp(-1, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,0,0) -@test Meshing.vertex_interp(1, Vec(0,0,0), Vec(0,1,0), -1, 1) == Vec(0,1,0) + @testset "marching tetrahedra" begin + m1 = @test_warn "is deprecated" GLNormalMesh(sdf) + m2 = @test_warn "is deprecated" GLNormalMesh(sdf, 1e-3) + m3 = @test_nowarn GLNormalMesh(sdf, MarchingTetrahedra()) + @test vertices(m1) == vertices(m2) == vertices(m3) + @test faces(m1) == faces(m2) == faces(m3) -# marching cubes -sdf = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v - sqrt(sum(dot(v,v))) - 1 # sphere + m4 = @test_warn "is deprecated" GLNormalMesh(sdf.data, 0.5) + m5 = @test_nowarn GLNormalMesh(sdf.data, MarchingTetrahedra(0.5)) + @test vertices(m4) == vertices(m5) + @test faces(m4) == faces(m5) + end + end end -m = marching_cubes(sdf,0) -m2 = marching_cubes(sdf) -@test length(vertices(m)) == 10968 -@test length(faces(m)) == 3656 -@test m == m2 if "--profile" in ARGS HomogenousMesh(s2)