diff --git a/src/algorithmtypes.jl b/src/algorithmtypes.jl index d21a58b..2f63c84 100644 --- a/src/algorithmtypes.jl +++ b/src/algorithmtypes.jl @@ -62,7 +62,7 @@ making this algorithm useful for topological analysis. - `iso` specifies the iso level to use for surface extraction. - `eps` is the tolerence around a voxel corner to ensure manifold mesh generation. -- `reduceverts` reserved for future use. +- `reduceverts` (default: true) if false, vertices will not be unique and have repeated face references. - `insidepositive` (default: false) set true if the sign convention inside the surface is positive, common for NRRD and DICOM data """ struct MarchingTetrahedra{T} <: AbstractMeshingAlgorithm diff --git a/src/marching_cubes.jl b/src/marching_cubes.jl index a34a00b..3a7322c 100644 --- a/src/marching_cubes.jl +++ b/src/marching_cubes.jl @@ -2,13 +2,6 @@ #Look up Table include("lut/mc.jl") -""" -Construct a `HomogenousMesh` from a `SignedDistanceField` using the -marching cubes algorithm. This method is faster than Marching Tetrahedra -and generates fewer vertices and faces (about 1/4 as many). -However it may generate non-manifold meshes, while Marching -Tetrahedra guarentees a manifold mesh. -""" function isosurface(sdf::AbstractArray{T, 3}, method::MarchingCubes, ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{3, Int}; origin=VertType(-1,-1,-1), widths=VertType(2,2,2)) where {T, VertType, FaceType} nx, ny, nz = size(sdf) diff --git a/src/marching_tetrahedra.jl b/src/marching_tetrahedra.jl index ae7bf6c..63de298 100644 --- a/src/marching_tetrahedra.jl +++ b/src/marching_tetrahedra.jl @@ -68,16 +68,25 @@ present. scale, origin, vts::Dict, vtsAry::Vector, - eps::Real) + eps::Real, + reduceverts::Bool) VertType = eltype(vtsAry) - vId = vertId(e, x, y, z, nx, ny) - - haskey(vts, vId) && return vts[vId] + if reduceverts + vId = vertId(e, x, y, z, nx, ny) + haskey(vts, vId) && return vts[vId] + end + # calculate vert position v = vertPos(e, x, y, z, scale, origin, vals, iso, eps, VertType) push!(vtsAry, v) - vts[vId] = length(vtsAry) + + # if deduplicting, push to dict + if reduceverts + vts[vId] = length(vtsAry) + end + + return length(vtsAry) end """ @@ -101,7 +110,7 @@ containers as necessary. """ function procVox(vals, iso::Real, x, y, z, nx, ny, scale, origin, vts::Dict, vtsAry::Vector, fcs::Vector, - eps::Real, cubeindex) + eps::Real, cubeindex, reduceverts) VertType = eltype(vtsAry) FaceType = eltype(fcs) # check each sub-tetrahedron in the voxel @@ -113,24 +122,19 @@ function procVox(vals, iso::Real, x, y, z, nx, ny, scale, origin, # add the face to the list push!(fcs, FaceType( - getVertId(voxEdgeId(i, e[1]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps), - getVertId(voxEdgeId(i, e[2]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps), - getVertId(voxEdgeId(i, e[3]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps))) + getVertId(voxEdgeId(i, e[1]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts), + getVertId(voxEdgeId(i, e[2]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts), + getVertId(voxEdgeId(i, e[3]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts))) # bail if there are no more faces iszero(e[4]) && continue push!(fcs, FaceType( - getVertId(voxEdgeId(i, e[4]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps), - getVertId(voxEdgeId(i, e[5]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps), - getVertId(voxEdgeId(i, e[6]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps))) + getVertId(voxEdgeId(i, e[4]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts), + getVertId(voxEdgeId(i, e[5]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts), + getVertId(voxEdgeId(i, e[6]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts))) end end - -""" -Given a 3D array and an isovalue, extracts a mesh represention of the -an approximate isosurface by the method of marching tetrahedra. -""" function isosurface(sdf::AbstractArray{T, 3}, method::MarchingTetrahedra, ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{3, Int}; origin=VertType(-1,-1,-1), widths=VertType(2,2,2)) where {T, VertType, FaceType} @@ -159,7 +163,7 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingTetrahedra, ::Type _no_triangles(cubeindex) && continue - procVox(vals, method.iso, i, j, k, nx, ny, scale, origin, vts, vtsAry, fcs, method.eps, cubeindex) + procVox(vals, method.iso, i, j, k, nx, ny, scale, origin, vts, vtsAry, fcs, method.eps, cubeindex, method.reduceverts) end vtsAry,fcs @@ -215,7 +219,7 @@ function isosurface(f::Function, method::MarchingTetrahedra, _no_triangles(cubeindex) && continue - procVox(vals, method.iso, i, j, k, nx, ny, scale, origin, vts, vtsAry, fcs, method.eps, cubeindex) + procVox(vals, method.iso, i, j, k, nx, ny, scale, origin, vts, vtsAry, fcs, method.eps, cubeindex, method.reduceverts) end vtsAry,fcs diff --git a/src/surface_nets.jl b/src/surface_nets.jl index 6b52fcc..4e6ed40 100644 --- a/src/surface_nets.jl +++ b/src/surface_nets.jl @@ -11,10 +11,6 @@ #Look up Table include("lut/sn.jl") -""" -Generate a mesh using naive surface nets. -This takes the center of mass of the voxel as the vertex for each cube. -""" function isosurface(sdf::AbstractArray{T, 3}, method::NaiveSurfaceNets, ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{4, Int}; origin=VertType(-1,-1,-1), widths=VertType(2,2,2)) where {T, VertType, FaceType} @@ -90,10 +86,6 @@ function isosurface(sdf::AbstractArray{T, 3}, method::NaiveSurfaceNets, ::Type{V vertices, faces end -""" -Generate a mesh using naive surface nets. -This takes the center of mass of the voxel as the vertex for each cube. -""" function isosurface(f::Function, method::NaiveSurfaceNets, ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{4, Int}; origin=VertType(-1,-1,-1), widths=VertType(2,2,2), diff --git a/test/runtests.jl b/test/runtests.jl index de7747f..6ec2aac 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -143,6 +143,21 @@ end @test length(faces(m)) == length(faces(mf)) end +@testset "marching tetrahedra" begin + a1 = MarchingTetrahedra() + a2 = MarchingTetrahedra(reduceverts=false) + m1 = SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), a1, samples=(21,21,21)) do v + sqrt(sum(dot(v,v))) - 1 # sphere + end + m2 = SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), a2, samples=(21,21,21)) do v + sqrt(sum(dot(v,v))) - 1 # sphere + end + @test length(m1.vertices) == 5582 + @test length(m2.vertices) == 33480 + @test length(m1.faces) == 11160 + @test length(m2.faces) == length(m1.faces) +end + @testset "sign flips" begin for algo_type in algos algo = algo_type() @@ -187,7 +202,10 @@ end resolution = 0.1 sdf = SignedDistanceField(f, bounds, resolution) - for algorithm in (MarchingCubes(0.5), MarchingTetrahedra(0.5)) + for algorithm in (MarchingCubes(0.5), + MarchingTetrahedra(0.5), + MarchingCubes(iso=0.5,reduceverts=false), + MarchingTetrahedra(iso=0.5,reduceverts=false)) @testset "$(typeof(algorithm))" begin mesh = @inferred GLNormalMesh(sdf, algorithm) # should be centered on the origin