Sjk/sn updates2

bench/benchmark.jl

@@ -15,7 +15,7 @@ println("Benchmarking Meshing.jl...")
 #
 
 algos_sdf = [MarchingCubes(), MarchingTetrahedra(), NaiveSurfaceNets()]
-algos_fn = [MarchingCubes(), NaiveSurfaceNets()]
+algos_fn = [MarchingCubes(), MarchingTetrahedra(), NaiveSurfaceNets()]
 
 #
 # Benchmark SDF constructon for SDF and Direct sample comparisons

docs/src/api.md

@@ -33,9 +33,10 @@ Three meshing algorithms exist:
 * `MarchingTetrahedra()`
 * `NaiveSurfaceNets()`
 
-Each takes an optional `iso` and `eps` parameter, e.g. `MarchingCubes(0.0,1e-6)`.
+Each takes optional `iso`, `eps`, and `insidepositive` parameters, e.g. `MarchingCubes(iso=0.0,eps=1e-6,insidepositive=false)`.
 
 Here `iso` controls the offset for the boundary detection. By default this is set to 0. `eps` is the detection tolerance for a voxel edge intersection.
+`insidepositive` sets the sign convention for inside/outside the surface (default: false).
 
 Users must construct an algorithm type and use it as an argument to a GeometryTypes mesh call or `isosurface` call.
 

src/algorithmtypes.jl

@@ -63,7 +63,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.
-- `insidepositive` reserved for future use.
+- `insidepositive` (default: false) set true if the sign convention inside the surface is positive, common for NRRD and DICOM data
 """
 struct MarchingTetrahedra{T} <: AbstractMeshingAlgorithm
     iso::T
@@ -86,7 +86,7 @@ however the algorithm does not guarantee accuracy and generates quad faces.
 - `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.
-- `insidepositive` reserved for future use.
+- `insidepositive` (default: false) set true if the sign convention inside the surface is positive, common for NRRD and DICOM data
 """
 struct NaiveSurfaceNets{T} <: AbstractMeshingAlgorithm
     iso::T

src/common.jl

@@ -39,6 +39,15 @@ where the sign convention indicates positive inside the surface
     cubeindex
 end
 
+"""
+    no_triangles(cubeindex)
+
+Called after `_get_cubeindex`. Determines if a voxel index has triangles.
+"""
+@inline function _no_triangles(cubeindex::UInt8)
+    cubeindex == 0x00 || cubeindex == 0xff
+end
+
 """
     _determine_types(meshtype, fieldtype=Float64, facelen=3)
 

src/lut/mt.jl

@@ -79,6 +79,15 @@ const subTets = ((0x01, 0x03, 0x02, 0x07),
                  (0x01, 0x02, 0x06, 0x07),
                  (0x01, 0x05, 0x08, 0x07),
                  (0x01, 0x06, 0x05, 0x07))
+
+# voxel corners that comprise each of the six tetrahedra, masking cubeindex
+const subTetsMask = ((0x04, 0x02),
+                     (0x80, 0x08),
+                     (0x08, 0x04),
+                     (0x02, 0x20),
+                     (0x10, 0x80),
+                     (0x20, 0x10))
+
 # tetrahedron corners for each edge (indices 1-4)
 const tetEdgeCrnrs = ((0x01, 0x02),
                       (0x02, 0x03),

src/lut/sn.jl

@@ -1,29 +1,35 @@
-const cube_edges = ( 0, 1, 0, 2, 0, 4, 1, 3, 1, 5, 2, 3,
-                     2, 6, 3, 7, 4, 5, 4, 6, 5, 7, 6, 7 )
-const sn_edge_table = (7, 25, 30, 98, 101, 123, 124, 168, 175, 177, 182, 202,
-                       205, 211, 212, 772, 771, 797, 794, 870, 865, 895, 888,
-                       940, 939, 949, 946, 974, 969, 983, 976, 1296, 1303, 1289,
-                       1294, 1394, 1397, 1387, 1388, 1464, 1471, 1441, 1446,
-                       1498, 1501, 1475, 1476, 1556, 1555, 1549, 1546, 1654,
-                       1649, 1647, 1640, 1724, 1723, 1701, 1698, 1758, 1753,
-                       1735, 1728, 2624, 2631, 2649, 2654, 2594, 2597, 2619,
-                       2620, 2792, 2799, 2801, 2806, 2698, 2701, 2707, 2708,
-                       2372, 2371, 2397, 2394, 2342, 2337, 2367, 2360, 2540,
-                       2539, 2549, 2546, 2446, 2441, 2455, 2448, 3920, 3927,
-                       3913, 3918, 3890, 3893, 3883, 3884, 4088, 4095, 4065,
-                       4070, 3994, 3997, 3971, 3972, 3156, 3155, 3149, 3146,
-                       3126, 3121, 3119, 3112, 3324, 3323, 3301, 3298, 3230,
-                       3225, 3207, 3200, 3200, 3207, 3225, 3230, 3298, 3301,
-                       3323, 3324, 3112, 3119, 3121, 3126, 3146, 3149, 3155,
-                       3156, 3972, 3971, 3997, 3994, 4070, 4065, 4095, 4088,
-                       3884, 3883, 3893, 3890, 3918, 3913, 3927, 3920, 2448,
-                       2455, 2441, 2446, 2546, 2549, 2539, 2540, 2360, 2367,
-                       2337, 2342, 2394, 2397, 2371, 2372, 2708, 2707, 2701,
-                       2698, 2806, 2801, 2799, 2792, 2620, 2619, 2597, 2594,
-                       2654, 2649, 2631, 2624, 1728, 1735, 1753, 1758, 1698,
-                       1701, 1723, 1724, 1640, 1647, 1649, 1654, 1546, 1549,
-                       1555, 1556, 1476, 1475, 1501, 1498, 1446, 1441, 1471,
-                       1464, 1388, 1387, 1397, 1394, 1294, 1289, 1303, 1296,
-                       976, 983, 969, 974, 946, 949, 939, 940, 888, 895, 865,
-                       870, 794, 797, 771, 772, 212, 211, 205, 202, 182, 177,
-                       175, 168, 124, 123, 101, 98, 30, 25, 7)
+const cube_edges = (0x00, 0x01, 0x00, 0x02, 0x00, 0x04, 0x01, 0x03, 0x01, 0x05, 0x02, 0x03,
+                    0x02, 0x06, 0x03, 0x07, 0x04, 0x05, 0x04, 0x06, 0x05, 0x07, 0x06, 0x07)
+
+const sn_edge_table = (0x0007, 0x0019, 0x001e, 0x0062, 0x0065, 0x007b, 0x007c, 0x00a8,
+                       0x00af, 0x00b1, 0x00b6, 0x00ca, 0x00cd, 0x00d3, 0x00d4, 0x0304,
+                       0x0303, 0x031d, 0x031a, 0x0366, 0x0361, 0x037f, 0x0378, 0x03ac,
+                       0x03ab, 0x03b5, 0x03b2, 0x03ce, 0x03c9, 0x03d7, 0x03d0, 0x0510,
+                       0x0517, 0x0509, 0x050e, 0x0572, 0x0575, 0x056b, 0x056c, 0x05b8,
+                       0x05bf, 0x05a1, 0x05a6, 0x05da, 0x05dd, 0x05c3, 0x05c4, 0x0614,
+                       0x0613, 0x060d, 0x060a, 0x0676, 0x0671, 0x066f, 0x0668, 0x06bc,
+                       0x06bb, 0x06a5, 0x06a2, 0x06de, 0x06d9, 0x06c7, 0x06c0, 0x0a40,
+                       0x0a47, 0x0a59, 0x0a5e, 0x0a22, 0x0a25, 0x0a3b, 0x0a3c, 0x0ae8,
+                       0x0aef, 0x0af1, 0x0af6, 0x0a8a, 0x0a8d, 0x0a93, 0x0a94, 0x0944,
+                       0x0943, 0x095d, 0x095a, 0x0926, 0x0921, 0x093f, 0x0938, 0x09ec,
+                       0x09eb, 0x09f5, 0x09f2, 0x098e, 0x0989, 0x0997, 0x0990, 0x0f50,
+                       0x0f57, 0x0f49, 0x0f4e, 0x0f32, 0x0f35, 0x0f2b, 0x0f2c, 0x0ff8,
+                       0x0fff, 0x0fe1, 0x0fe6, 0x0f9a, 0x0f9d, 0x0f83, 0x0f84, 0x0c54,
+                       0x0c53, 0x0c4d, 0x0c4a, 0x0c36, 0x0c31, 0x0c2f, 0x0c28, 0x0cfc,
+                       0x0cfb, 0x0ce5, 0x0ce2, 0x0c9e, 0x0c99, 0x0c87, 0x0c80, 0x0c80,
+                       0x0c87, 0x0c99, 0x0c9e, 0x0ce2, 0x0ce5, 0x0cfb, 0x0cfc, 0x0c28,
+                       0x0c2f, 0x0c31, 0x0c36, 0x0c4a, 0x0c4d, 0x0c53, 0x0c54, 0x0f84,
+                       0x0f83, 0x0f9d, 0x0f9a, 0x0fe6, 0x0fe1, 0x0fff, 0x0ff8, 0x0f2c,
+                       0x0f2b, 0x0f35, 0x0f32, 0x0f4e, 0x0f49, 0x0f57, 0x0f50, 0x0990,
+                       0x0997, 0x0989, 0x098e, 0x09f2, 0x09f5, 0x09eb, 0x09ec, 0x0938,
+                       0x093f, 0x0921, 0x0926, 0x095a, 0x095d, 0x0943, 0x0944, 0x0a94,
+                       0x0a93, 0x0a8d, 0x0a8a, 0x0af6, 0x0af1, 0x0aef, 0x0ae8, 0x0a3c,
+                       0x0a3b, 0x0a25, 0x0a22, 0x0a5e, 0x0a59, 0x0a47, 0x0a40, 0x06c0,
+                       0x06c7, 0x06d9, 0x06de, 0x06a2, 0x06a5, 0x06bb, 0x06bc, 0x0668,
+                       0x066f, 0x0671, 0x0676, 0x060a, 0x060d, 0x0613, 0x0614, 0x05c4,
+                       0x05c3, 0x05dd, 0x05da, 0x05a6, 0x05a1, 0x05bf, 0x05b8, 0x056c,
+                       0x056b, 0x0575, 0x0572, 0x050e, 0x0509, 0x0517, 0x0510, 0x03d0,
+                       0x03d7, 0x03c9, 0x03ce, 0x03b2, 0x03b5, 0x03ab, 0x03ac, 0x0378,
+                       0x037f, 0x0361, 0x0366, 0x031a, 0x031d, 0x0303, 0x0304, 0x00d4,
+                       0x00d3, 0x00cd, 0x00ca, 0x00b6, 0x00b1, 0x00af, 0x00a8, 0x007c,
+                       0x007b, 0x0065, 0x0062, 0x001e, 0x0019, 0x0007)

src/marching_cubes.jl

@@ -41,7 +41,7 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingCubes, ::Type{Vert
         cubeindex = method.insidepositive ? _get_cubeindex_pos(iso_vals, method.iso) : _get_cubeindex(iso_vals, method.iso)
 
         # Cube is entirely in/out of the surface
-        (cubeindex == 0x00 || cubeindex == 0xff) && continue
+        _no_triangles(cubeindex) && continue
 
         points = mc_vert_points(xi,yi,zi,s,origin,VertType)
 
@@ -72,33 +72,38 @@ function isosurface(f::Function, method::MarchingCubes, ::Type{VertType}=SVector
     method.reduceverts && sizehint!(vts, mt*mt*5)
     !method.reduceverts && sizehint!(vts, mt*mt*6)
     sizehint!(fcs, mt*mt*2)
-    iso_vals = Vector{eltype(VertType)}(undef,8)
+    zv = zero(eltype(VertType))
+    iso_vals = (zv,zv,zv,zv,zv,zv,zv,zv)
     @inbounds for xi = 1:nx-1, yi = 1:ny-1, zi = 1:nz-1
 
-
         points = mc_vert_points(xi,yi,zi,s,origin,VertType)
 
         if zi == 1
-            for i = 1:8
-                iso_vals[i] = f(points[i])
-            end
+            iso_vals = (f(points[1]),
+                        f(points[2]),
+                        f(points[3]),
+                        f(points[4]),
+                        f(points[5]),
+                        f(points[6]),
+                        f(points[7]),
+                        f(points[8]))
         else
-            iso_vals[1] = iso_vals[5]
-            iso_vals[2] = iso_vals[6]
-            iso_vals[3] = iso_vals[7]
-            iso_vals[4] = iso_vals[8]
-            iso_vals[5] = f(points[5])
-            iso_vals[6] = f(points[6])
-            iso_vals[7] = f(points[7])
-            iso_vals[8] = f(points[8])
+            iso_vals = (iso_vals[5],
+                        iso_vals[6],
+                        iso_vals[7],
+                        iso_vals[8],
+                        f(points[5]),
+                        f(points[6]),
+                        f(points[7]),
+                        f(points[8]))
         end
 
         #Determine the index into the edge table which
         #tells us which vertices are inside of the surface
         cubeindex = method.insidepositive ? _get_cubeindex_pos(iso_vals, method.iso) : _get_cubeindex(iso_vals, method.iso)
 
         # Cube is entirely in/out of the surface
-        (cubeindex == 0x00 || cubeindex == 0xff) && continue
+        _no_triangles(cubeindex) && continue
 
         # Find the vertices where the surface intersects the cube
         # TODO this can use the underlying function to find the zero.
@@ -119,7 +124,7 @@ end
 
 Create triangles by adding every point within a triangle to the vertex vector.
 """
-@inline function _mc_create_triangles!(vts, fcs, vertlist, cubeindex, FaceType)
+@inline function _mc_create_triangles!(vts, fcs, vertlist, cubeindex, ::Type{FaceType}) where {FaceType}
     fct = length(vts) + 3
 
     push!(vts, vertlist[tri_table[cubeindex][1]],
@@ -162,7 +167,7 @@ end
 Create triangles by only adding unique vertices within the voxel.
 Each face may share a reference to a vertex with another face.
 """
-@inline function _mc_unique_triangles!(vts, fcs, vertlist, cubeindex, FaceType)
+@inline function _mc_unique_triangles!(vts, fcs, vertlist, cubeindex, ::Type{FaceType}) where {FaceType}
     fct = length(vts)
 
     vert_to_add = _mc_verts[cubeindex]
@@ -176,6 +181,7 @@ Each face may share a reference to a vertex with another face.
         elt = vert_to_add[i]
         push!(vts, vertlist[elt])
     end
+
     offsets = _mc_connectivity[cubeindex]
 
     # There is atleast one face so we can push it immediately
@@ -250,4 +256,4 @@ Returns a tuple of 8 points corresponding to each corner of a cube
      VertType(xi  ,yi-1,zi  ) .* scale .+ origin,
      VertType(xi  ,yi  ,zi  ) .* scale .+ origin,
      VertType(xi-1,yi  ,zi  ) .* scale .+ origin)
-end
+end