Merge pull request #60 from JuliaGeometry/sjk/mt_noreduce

[mt] add option to not hash vertices, useful for octtrees

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

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)

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

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),

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