From c99bcaf9a3608da3e9cde61f33db56915bb9ad12 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Fri, 20 Dec 2019 15:43:06 -0500 Subject: [PATCH 01/19] use tuples for isovals in function MC --- src/marching_cubes.jl | 32 +++++++++++++++++++------------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/src/marching_cubes.jl b/src/marching_cubes.jl index 654cb05..a6c8ee1 100644 --- a/src/marching_cubes.jl +++ b/src/marching_cubes.jl @@ -72,25 +72,31 @@ 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 @@ -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 \ No newline at end of file +end From 0495b3a6766f32496c04a29cccd08016b73367c9 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Fri, 20 Dec 2019 17:05:04 -0500 Subject: [PATCH 02/19] general cleanup to SN implementation --- src/common.jl | 9 +++++++++ src/surface_nets.jl | 24 ++++++++++-------------- 2 files changed, 19 insertions(+), 14 deletions(-) diff --git a/src/common.jl b/src/common.jl index 039822b..2bed6e0 100644 --- a/src/common.jl +++ b/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`. +""" +@inline function _no_triangles(cubeindex::UInt8) + cubeindex == 0x00 || cubeindex == 0xff +end + """ _determine_types(meshtype, fieldtype=Float64, facelen=3) diff --git a/src/surface_nets.jl b/src/surface_nets.jl index 104414a..3270f06 100644 --- a/src/surface_nets.jl +++ b/src/surface_nets.jl @@ -72,7 +72,7 @@ function isosurface(sdf::AbstractArray{T, 3}, method::NaiveSurfaceNets, ::Type{V mask = _get_cubeindex(grid, 0) # Check for early termination if cell does not intersect boundary - if mask == 0x00 || mask == 0xff + if _no_triangles(mask) xi += 1 n += 1 m += 1 @@ -82,7 +82,7 @@ function isosurface(sdf::AbstractArray{T, 3}, method::NaiveSurfaceNets, ::Type{V #Sum up edge intersections edge_mask = sn_edge_table[mask] - _sn_add_verts!(xi, yi, zi, vertices, grid, edge_mask, buffer, m, scale, origin, method.eps, Val(true), VertType) + _sn_add_verts!(xi, yi, zi, vertices, grid, edge_mask, buffer, m, scale, origin, method.eps, true, VertType) #Now we need to add faces together, to do this we just loop over 3 basis components x = (xi,yi,zi) @@ -210,7 +210,7 @@ function isosurface(f::Function, method::NaiveSurfaceNets, mask = _get_cubeindex(grid, 0) # Check for early termination if cell does not intersect boundary - if mask == 0x00 || mask == 0xff + if _no_triangles(mask) xi += 1 n += 1 m += 1 @@ -221,7 +221,7 @@ function isosurface(f::Function, method::NaiveSurfaceNets, edge_mask = sn_edge_table[mask] # add vertices - _sn_add_verts!(xi, yi, zi, vertices, grid, edge_mask, buffer, m, scale, origin, method.eps, Val(false), VertType) + _sn_add_verts!(xi, yi, zi, vertices, grid, edge_mask, buffer, m, scale, origin, method.eps, false, VertType) #Now we need to add faces together, to do this we just loop over 3 basis components x = (xi,yi,zi) @@ -278,9 +278,7 @@ end @inbounds for i=0:11 #Use edge mask to check if it is crossed - if (edge_mask & (1< eps - t = g0 / t - else - continue - end + + abs(t) <= eps && continue + + t = g0 / t #Interpolate vertices and add up intersections (this can be done without multiplying) # TODO lut table change may have made this incorrect @@ -319,7 +316,7 @@ end #Now we just average the edge intersections and add them to coordinate s = 1.0 / e_count - if typeof(translate_pt) === Val{true} + if translate_pt v = (VertType(xi,yi,zi) .+ s .* v) .* scale + origin else v = (VertType(xi,yi,zi) .+ s .* v)# * scale[i] + origin[i] @@ -329,4 +326,3 @@ end buffer[m+1] = length(vertices) push!(vertices, v) end - From d6eff5cabc3f4871f0ddd684c45fdf70cb966e0c Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Fri, 20 Dec 2019 20:34:30 -0500 Subject: [PATCH 03/19] cleanup vertex generation in SN --- src/lut/sn.jl | 5 +++-- src/surface_nets.jl | 33 +++++++++++++-------------------- 2 files changed, 16 insertions(+), 22 deletions(-) diff --git a/src/lut/sn.jl b/src/lut/sn.jl index 55f8b4c..93c22ee 100644 --- a/src/lut/sn.jl +++ b/src/lut/sn.jl @@ -1,5 +1,6 @@ -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 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 = (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, diff --git a/src/surface_nets.jl b/src/surface_nets.jl index 3270f06..5bc91c8 100644 --- a/src/surface_nets.jl +++ b/src/surface_nets.jl @@ -38,7 +38,7 @@ function isosurface(sdf::AbstractArray{T, 3}, method::NaiveSurfaceNets, ::Type{V sizehint!(faces,ceil(Int,maximum(dims)^2)) n = 0 - R = Array{Int}([1, (dims[1]+1), (dims[1]+1)*(dims[2]+1)]) + R = Array{Int}([1, (dims[1]+1), (dims[1]+1)*(dims[2]+1)]) #TODO buf_no = 1 buffer = fill(zero(Int),R[3]*2) @@ -275,7 +275,7 @@ end e_count = 0 #For every edge of the cube... - @inbounds for i=0:11 + @inbounds for i=0x00:0x0b #Use edge mask to check if it is crossed iszero(edge_mask & (1< Date: Fri, 20 Dec 2019 20:37:24 -0500 Subject: [PATCH 04/19] use _no_triangles in a few parts --- src/common.jl | 2 +- src/marching_cubes.jl | 6 +++--- src/surface_nets.jl | 20 +++++++++++--------- 3 files changed, 15 insertions(+), 13 deletions(-) diff --git a/src/common.jl b/src/common.jl index 2bed6e0..d1c8d62 100644 --- a/src/common.jl +++ b/src/common.jl @@ -42,7 +42,7 @@ end """ no_triangles(cubeindex) -Called after `_get_cubeindex`. +Called after `_get_cubeindex`. Determines if a voxel index has triangles. """ @inline function _no_triangles(cubeindex::UInt8) cubeindex == 0x00 || cubeindex == 0xff diff --git a/src/marching_cubes.jl b/src/marching_cubes.jl index a6c8ee1..2233c0c 100644 --- a/src/marching_cubes.jl +++ b/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) @@ -76,7 +76,6 @@ function isosurface(f::Function, method::MarchingCubes, ::Type{VertType}=SVector 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 @@ -104,7 +103,7 @@ function isosurface(f::Function, method::MarchingCubes, ::Type{VertType}=SVector 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. @@ -182,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 diff --git a/src/surface_nets.jl b/src/surface_nets.jl index 5bc91c8..bd09cbf 100644 --- a/src/surface_nets.jl +++ b/src/surface_nets.jl @@ -58,6 +58,8 @@ function isosurface(sdf::AbstractArray{T, 3}, method::NaiveSurfaceNets, ::Type{V xi=0 @inbounds while xi < dims[1]-1 + inds = (xi,yi,zi) + # Read in 8 field values around this vertex and store them in an array # Also calculate 8-bit mask, like in marching cubes, so we can speed up sign checks later @inbounds grid = (data[n+1], @@ -82,10 +84,9 @@ function isosurface(sdf::AbstractArray{T, 3}, method::NaiveSurfaceNets, ::Type{V #Sum up edge intersections edge_mask = sn_edge_table[mask] - _sn_add_verts!(xi, yi, zi, vertices, grid, edge_mask, buffer, m, scale, origin, method.eps, true, VertType) + _sn_add_verts!(inds, vertices, grid, edge_mask, buffer, m, scale, origin, method.eps, true, VertType) #Now we need to add faces together, to do this we just loop over 3 basis components - x = (xi,yi,zi) for i=0:2 #The first three entries of the edge_mask count the crossings along the edge if (edge_mask & (1< Date: Fri, 20 Dec 2019 21:59:43 -0500 Subject: [PATCH 05/19] [SN] factor out face generation --- src/surface_nets.jl | 106 +++++++++++++++++--------------------------- 1 file changed, 41 insertions(+), 65 deletions(-) diff --git a/src/surface_nets.jl b/src/surface_nets.jl index bd09cbf..a73f1c6 100644 --- a/src/surface_nets.jl +++ b/src/surface_nets.jl @@ -87,42 +87,8 @@ function isosurface(sdf::AbstractArray{T, 3}, method::NaiveSurfaceNets, ::Type{V _sn_add_verts!(inds, vertices, grid, edge_mask, buffer, m, scale, origin, method.eps, true, VertType) #Now we need to add faces together, to do this we just loop over 3 basis components - for i=0:2 - #The first three entries of the edge_mask count the crossings along the edge - if (edge_mask & (1< Date: Fri, 20 Dec 2019 22:02:11 -0500 Subject: [PATCH 06/19] [sn] use iszero more places --- src/surface_nets.jl | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/surface_nets.jl b/src/surface_nets.jl index a73f1c6..0fa5bb2 100644 --- a/src/surface_nets.jl +++ b/src/surface_nets.jl @@ -262,25 +262,21 @@ end function _sn_add_faces!(inds, faces, edge_mask, mask, buffer, m, R, ::Type{FaceType}) where {FaceType} for i=0:2 #The first three entries of the edge_mask count the crossings along the edge - if (edge_mask & (1< Date: Fri, 20 Dec 2019 21:20:12 -0500 Subject: [PATCH 07/19] [sn] use 16bit uints in lut --- src/lut/sn.jl | 59 ++++++++++++++++++++++++++++----------------------- 1 file changed, 32 insertions(+), 27 deletions(-) diff --git a/src/lut/sn.jl b/src/lut/sn.jl index 93c22ee..610c62f 100644 --- a/src/lut/sn.jl +++ b/src/lut/sn.jl @@ -1,30 +1,35 @@ 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 = (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 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) From 1c3b5b45ad977396ab80c750b2e435111e7f3754 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Fri, 20 Dec 2019 22:27:39 -0500 Subject: [PATCH 08/19] [sn] better masking logic --- src/surface_nets.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/surface_nets.jl b/src/surface_nets.jl index 0fa5bb2..5272999 100644 --- a/src/surface_nets.jl +++ b/src/surface_nets.jl @@ -217,7 +217,7 @@ end @inbounds for i=0x00:0x0b #Use edge mask to check if it is crossed - iszero(edge_mask & (1< Date: Fri, 20 Dec 2019 22:52:53 -0500 Subject: [PATCH 09/19] [sn] use ntuples in grid updates --- src/surface_nets.jl | 30 ++++++++++++++++++------------ 1 file changed, 18 insertions(+), 12 deletions(-) diff --git a/src/surface_nets.jl b/src/surface_nets.jl index 5272999..5524986 100644 --- a/src/surface_nets.jl +++ b/src/surface_nets.jl @@ -130,7 +130,8 @@ function isosurface(f::Function, method::NaiveSurfaceNets, buffer = fill(zero(Int),R[3]*2) - grid = Vector{eltype(VertType)}(undef,8) + zv = zero(eltype(VertType)) + grid = (zv,zv,zv,zv,zv,zv,zv,zv) #March over the voxel grid zi = 0 @@ -160,18 +161,23 @@ function isosurface(f::Function, method::NaiveSurfaceNets, VertType(xi+1,yi+1,zi+1).* scale + origin) if xi == 0 - for i = 1:8 - grid[i] = f(points[i]) - end + grid = (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 - grid[1] = grid[2] - grid[2] = f(points[2]) - grid[3] = grid[4] - grid[4] = f(points[4]) - grid[5] = grid[6] - grid[6] = f(points[6]) - grid[7] = grid[8] - grid[8] = f(points[8]) + grid = (grid[2], + f(points[2]), + grid[4], + f(points[4]), + grid[6], + f(points[6]), + grid[8], + f(points[8])) end # Also calculate 8-bit mask, like in marching cubes, so we can speed up sign checks later From 3b5ca0e48cf38411bf6e953bf2b2e57f2803bdb7 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Fri, 20 Dec 2019 23:12:02 -0500 Subject: [PATCH 10/19] [sn] use for loops where possible --- src/surface_nets.jl | 27 ++++++--------------------- 1 file changed, 6 insertions(+), 21 deletions(-) diff --git a/src/surface_nets.jl b/src/surface_nets.jl index 5524986..c36e895 100644 --- a/src/surface_nets.jl +++ b/src/surface_nets.jl @@ -52,11 +52,8 @@ function isosurface(sdf::AbstractArray{T, 3}, method::NaiveSurfaceNets, ::Type{V # The contents of the buffer will be the indices of the vertices on the previous x/y slice of the volume m = 1 + (dims[1]+1) * (1 + buf_no * (dims[2]+1)) - yi=0 - @inbounds while yi Date: Sat, 21 Dec 2019 00:00:49 -0500 Subject: [PATCH 11/19] [sn] improve point creation --- src/surface_nets.jl | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/src/surface_nets.jl b/src/surface_nets.jl index c36e895..f6cd355 100644 --- a/src/surface_nets.jl +++ b/src/surface_nets.jl @@ -139,15 +139,15 @@ function isosurface(f::Function, method::NaiveSurfaceNets, inds = (xi,yi,zi) - # Read in 8 field values around this vertex and store them in an array - points = (VertType(xi,yi,zi).* scale + origin, - VertType(xi+1,yi,zi).* scale + origin, - VertType(xi,yi+1,zi).* scale + origin, - VertType(xi+1,yi+1,zi).* scale + origin, - VertType(xi,yi,zi+1).* scale + origin, - VertType(xi+1,yi,zi+1).* scale + origin, - VertType(xi,yi+1,zi+1).* scale + origin, - VertType(xi+1,yi+1,zi+1).* scale + origin) + if xi == 0 + points = (VertType(xi,yi,zi).* scale + origin, + VertType(xi+1,yi,zi).* scale + origin, + VertType(xi,yi+1,zi).* scale + origin, + VertType(xi+1,yi+1,zi).* scale + origin, + VertType(xi,yi,zi+1).* scale + origin, + VertType(xi+1,yi,zi+1).* scale + origin, + VertType(xi,yi+1,zi+1).* scale + origin, + VertType(xi+1,yi+1,zi+1).* scale + origin) if xi == 0 grid = (f(points[1]), @@ -159,6 +159,14 @@ function isosurface(f::Function, method::NaiveSurfaceNets, f(points[7]), f(points[8])) else + points = (points[2], + VertType(xi+1,yi,zi).* scale + origin, + points[4], + VertType(xi+1,yi+1,zi).* scale + origin, + points[6], + VertType(xi+1,yi,zi+1).* scale + origin, + points[8], + VertType(xi+1,yi+1,zi+1).* scale + origin) grid = (grid[2], f(points[2]), grid[4], From 0d22748379016c6d29b118d18cfc55fc098e4795 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Sat, 21 Dec 2019 00:01:41 -0500 Subject: [PATCH 12/19] [sn] improve sdf indexing and isolevel handling --- src/surface_nets.jl | 45 +++++++++++++++++++-------------------------- 1 file changed, 19 insertions(+), 26 deletions(-) diff --git a/src/surface_nets.jl b/src/surface_nets.jl index f6cd355..7a422b0 100644 --- a/src/surface_nets.jl +++ b/src/surface_nets.jl @@ -21,15 +21,6 @@ function isosurface(sdf::AbstractArray{T, 3}, method::NaiveSurfaceNets, ::Type{V scale = widths ./ VertType(size(sdf) .- 1) # subtract 1 because an SDF with N points per side has N-1 cells dims = size(sdf) - data = vec(sdf) - - # Run iso surface additions here - # TODO - if method.iso != 0.0 - for i = eachindex(data) - data[i] -= method.iso - end - end vertices = VertType[] faces = FaceType[] @@ -38,7 +29,7 @@ function isosurface(sdf::AbstractArray{T, 3}, method::NaiveSurfaceNets, ::Type{V sizehint!(faces,ceil(Int,maximum(dims)^2)) n = 0 - R = Array{Int}([1, (dims[1]+1), (dims[1]+1)*(dims[2]+1)]) #TODO + R = Array{Int}([1, (dims[1]+1), (dims[1]+1)*(dims[2]+1)]) buf_no = 1 buffer = fill(zero(Int),R[3]*2) @@ -59,24 +50,26 @@ function isosurface(sdf::AbstractArray{T, 3}, method::NaiveSurfaceNets, ::Type{V # Read in 8 field values around this vertex and store them in an array # Also calculate 8-bit mask, like in marching cubes, so we can speed up sign checks later - @inbounds grid = (data[n+1], - data[n+2], - data[n+dims[1]+1], - data[n+dims[1]+2], - data[n+dims[1]*2+1 + dims[1]*(dims[2]-2)], - data[n+dims[1]*2+2 + dims[1]*(dims[2]-2)], - data[n+dims[1]*3+1 + dims[1]*(dims[2]-2)], - data[n+dims[1]*3+2 + dims[1]*(dims[2]-2)]) + @inbounds grid = (sdf[xi+1,yi+1,zi+1], + sdf[xi+2,yi+1,zi+1], + sdf[xi+1,yi+2,zi+1], + sdf[xi+2,yi+2,zi+1], + sdf[xi+1,yi+1,zi+2], + sdf[xi+2,yi+1,zi+2], + sdf[xi+1,yi+2,zi+2], + sdf[xi+2,yi+2,zi+2]) - mask = _get_cubeindex(grid, 0) + mask = _get_cubeindex(grid, method.iso) # Check for early termination if cell does not intersect boundary if _no_triangles(mask) - n += 1 m += 1 continue end + # iso level correction + grid = grid .- method.iso + #Sum up edge intersections edge_mask = sn_edge_table[mask] @@ -85,14 +78,11 @@ function isosurface(sdf::AbstractArray{T, 3}, method::NaiveSurfaceNets, ::Type{V #Now we need to add faces together, to do this we just loop over 3 basis components _sn_add_faces!(inds, faces, edge_mask, mask, buffer, m, R, FaceType) - n += 1 m += 1 end - n += 1 m += 2 end zi += 1 - n+=dims[1] buf_no = xor(buf_no,1) R[3]=-R[3] end @@ -123,7 +113,9 @@ function isosurface(f::Function, method::NaiveSurfaceNets, buffer = fill(zero(Int),R[3]*2) zv = zero(eltype(VertType)) + zvt = zero(VertType) grid = (zv,zv,zv,zv,zv,zv,zv,zv) + points = (zvt,zvt,zvt,zvt,zvt,zvt,zvt,zvt) #March over the voxel grid zi = 0 @@ -148,8 +140,6 @@ function isosurface(f::Function, method::NaiveSurfaceNets, VertType(xi+1,yi,zi+1).* scale + origin, VertType(xi,yi+1,zi+1).* scale + origin, VertType(xi+1,yi+1,zi+1).* scale + origin) - - if xi == 0 grid = (f(points[1]), f(points[2]), f(points[3]), @@ -178,7 +168,7 @@ function isosurface(f::Function, method::NaiveSurfaceNets, end # Also calculate 8-bit mask, like in marching cubes, so we can speed up sign checks later - mask = _get_cubeindex(grid, 0) + mask = _get_cubeindex(grid, method.iso) # Check for early termination if cell does not intersect boundary if _no_triangles(mask) @@ -186,6 +176,9 @@ function isosurface(f::Function, method::NaiveSurfaceNets, continue end + # iso level correction + grid = grid .- method.iso + #Sum up edge intersections edge_mask = sn_edge_table[mask] From d27e0063884cf4b89f49db679ad32548b044a6e2 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Sat, 21 Dec 2019 00:13:14 -0500 Subject: [PATCH 13/19] [mc] fix type signatures on triangle creations --- src/marching_cubes.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/marching_cubes.jl b/src/marching_cubes.jl index 2233c0c..abce8d8 100644 --- a/src/marching_cubes.jl +++ b/src/marching_cubes.jl @@ -124,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]], @@ -167,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] From 6b3c20262e6b66a0f13768c0dba720543049a592 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Sat, 21 Dec 2019 20:10:24 -0500 Subject: [PATCH 14/19] [mt] use cubeindex and new map in vertID to help support for method.insidepositive --- src/lut/mt.jl | 9 +++++++++ src/marching_tetrahedra.jl | 24 ++++++++++-------------- 2 files changed, 19 insertions(+), 14 deletions(-) diff --git a/src/lut/mt.jl b/src/lut/mt.jl index f787da7..ae0e280 100644 --- a/src/lut/mt.jl +++ b/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), diff --git a/src/marching_tetrahedra.jl b/src/marching_tetrahedra.jl index b1848b1..0d33816 100644 --- a/src/marching_tetrahedra.jl +++ b/src/marching_tetrahedra.jl @@ -4,19 +4,16 @@ include("lut/mt.jl") """ - tetIx(tIx, vals, iso::Real) + tetIx(tIx, cubeindex) Determines which case in the triangle table we are dealing with """ -@inline function tetIx(tIx, vals, iso::Real) - v1 = vals[subTets[tIx][1]] - v2 = vals[subTets[tIx][2]] - v3 = vals[subTets[tIx][3]] - v4 = vals[subTets[tIx][4]] - ix = v1 < iso ? 0x01 : 0x00 - (v2 < iso) && (ix |= 0x02) - (v3 < iso) && (ix |= 0x04) - (v4 < iso) && (ix |= 0x08) +@inline function tetIx(tIx, cubeindex) + @inbounds v1 = subTetsMask[tIx][1] + @inbounds v2 = subTetsMask[tIx][2] + ix = 0x01 & cubeindex | (0x40 & cubeindex) >> 0x03 + !iszero(v1 & cubeindex) && (ix |= 0x02) + !iszero(v2 & cubeindex) && (ix |= 0x04) ix end @@ -103,12 +100,12 @@ containers as necessary. """ function procVox(vals, iso::Real, x, y, z, nx, ny, vts::Dict, vtsAry::Vector, fcs::Vector, - eps::Real) + eps::Real, cubeindex) VertType = eltype(vtsAry) FaceType = eltype(fcs) # check each sub-tetrahedron in the voxel @inbounds for i = 1:6 - tIx = tetIx(i, vals, iso) + tIx = tetIx(i, cubeindex) (tIx == 0x00 || tIx == 0x0f) && continue e = tetTri[tIx] @@ -154,7 +151,7 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingTetrahedra, ::Type sdf[i+1, j, k+1]) cubeindex = _get_cubeindex(vals,method.iso) if cubeindex != 0x00 && cubeindex != 0xff - procVox(vals, method.iso, i, j, k, nx, ny, vts, vtsAry, fcs, method.eps) + procVox(vals, method.iso, i, j, k, nx, ny, vts, vtsAry, fcs, method.eps, cubeindex) end end @@ -176,4 +173,3 @@ function _correct_vertices!(vts, size, origin, widths, ::Type{VertType}) where { vts[i] = (vts[i] .- 1) .* s .+ origin # subtract 1 to fix 1-indexing end end - From 5cbc64a6090fb3a8ace21087314edf71ea026d3b Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Sat, 21 Dec 2019 22:18:58 -0500 Subject: [PATCH 15/19] [mt] run correct vertices in array push --- src/marching_tetrahedra.jl | 41 +++++++++++++------------------------- 1 file changed, 14 insertions(+), 27 deletions(-) diff --git a/src/marching_tetrahedra.jl b/src/marching_tetrahedra.jl index 0d33816..c1c1bca 100644 --- a/src/marching_tetrahedra.jl +++ b/src/marching_tetrahedra.jl @@ -39,7 +39,8 @@ occurs. eps represents the "bump" factor to keep vertices away from voxel corners (thereby preventing degeneracies). """ -@inline function vertPos(e, x, y, z, vals::NTuple{8,T}, iso::Real, eps::Real, ::Type{VertType}) where {T<:Real, VertType} +@inline function vertPos(e, x, y, z, scale, origin, vals::V, iso, eps, ::Type{VertType}) where {V, VertType} + T = eltype(vals) ixs = voxEdgeCrnrs[e] srcVal = vals[ixs[1]] @@ -49,7 +50,7 @@ corners (thereby preventing degeneracies). c1 = voxCrnrPos(VertType)[ixs[1]] c2 = voxCrnrPos(VertType)[ixs[2]] - VertType(x,y,z) + c1 .* b + c2.* a + ((VertType(x,y,z) + c1 .* b + c2.* a) .- 1) .* scale .+ origin end """ @@ -64,6 +65,7 @@ present. """ @inline function getVertId(e, x, y, z, nx, ny, vals, iso::Real, + scale, origin, vts::Dict, vtsAry::Vector, eps::Real) @@ -72,7 +74,7 @@ present. vId = vertId(e, x, y, z, nx, ny) # TODO we can probably immediately construct the vertex array here and use vert id to map to sequential ordering if !haskey(vts, vId) - v = vertPos(e, x, y, z, vals, iso, eps, VertType) + v = vertPos(e, x, y, z, scale, origin, vals, iso, eps, VertType) push!(vtsAry, v) vts[vId] = length(vtsAry) end @@ -98,7 +100,7 @@ end Processes a voxel, adding any new vertices and faces to the given containers as necessary. """ -function procVox(vals, iso::Real, x, y, z, nx, ny, +function procVox(vals, iso::Real, x, y, z, nx, ny, scale, origin, vts::Dict, vtsAry::Vector, fcs::Vector, eps::Real, cubeindex) VertType = eltype(vtsAry) @@ -112,16 +114,16 @@ function procVox(vals, iso::Real, x, y, z, nx, ny, # add the face to the list push!(fcs, FaceType( - getVertId(voxEdgeId(i, e[1]), x, y, z, nx, ny, vals, iso, vts, vtsAry, eps), - getVertId(voxEdgeId(i, e[2]), x, y, z, nx, ny, vals, iso, vts, vtsAry, eps), - getVertId(voxEdgeId(i, e[3]), x, y, z, nx, ny, vals, iso, vts, vtsAry, eps))) + 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))) # 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, vts, vtsAry, eps), - getVertId(voxEdgeId(i, e[5]), x, y, z, nx, ny, vals, iso, vts, vtsAry, eps), - getVertId(voxEdgeId(i, e[6]), x, y, z, nx, ny, vals, iso, vts, vtsAry, eps))) + 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))) end end @@ -139,6 +141,7 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingTetrahedra, ::Type sizehint!(vtsAry, div(length(sdf),8)) sizehint!(fcs, div(length(sdf),4)) # process each voxel + scale = widths ./ VertType(size(sdf) .- 1) nx::Int, ny::Int, nz::Int = size(sdf) @inbounds for k = 1:nz-1, j = 1:ny-1, i= 1:nx-1 vals = (sdf[i, j, k], @@ -151,25 +154,9 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingTetrahedra, ::Type sdf[i+1, j, k+1]) cubeindex = _get_cubeindex(vals,method.iso) if cubeindex != 0x00 && cubeindex != 0xff - procVox(vals, method.iso, i, j, k, nx, ny, vts, vtsAry, fcs, method.eps, cubeindex) + procVox(vals, method.iso, i, j, k, nx, ny, scale, origin, vts, vtsAry, fcs, method.eps, cubeindex) end end - _correct_vertices!(vtsAry, size(sdf), origin, widths, VertType) - vtsAry,fcs end - -""" - _correct_vertices!(vts, size, origin, widths, VertType) - -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, size, origin, widths, ::Type{VertType}) where {VertType} - s = widths ./ VertType(size .- 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 .+ origin # subtract 1 to fix 1-indexing - end -end From a15fcb985a80512df64a65fb4ea506207c571c6a Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Sat, 21 Dec 2019 22:38:12 -0500 Subject: [PATCH 16/19] [mt] formatting updates --- src/marching_tetrahedra.jl | 46 +++++++++++++++++++++----------------- 1 file changed, 25 insertions(+), 21 deletions(-) diff --git a/src/marching_tetrahedra.jl b/src/marching_tetrahedra.jl index c1c1bca..907f76f 100644 --- a/src/marching_tetrahedra.jl +++ b/src/marching_tetrahedra.jl @@ -72,13 +72,12 @@ present. VertType = eltype(vtsAry) vId = vertId(e, x, y, z, nx, ny) - # TODO we can probably immediately construct the vertex array here and use vert id to map to sequential ordering - if !haskey(vts, vId) - v = vertPos(e, x, y, z, scale, origin, vals, iso, eps, VertType) - push!(vtsAry, v) - vts[vId] = length(vtsAry) - end - vts[vId] + + haskey(vts, vId) && return vts[vId] + + v = vertPos(e, x, y, z, scale, origin, vals, iso, eps, VertType) + push!(vtsAry, v) + vts[vId] = length(vtsAry) end """ @@ -134,8 +133,9 @@ an approximate isosurface by the method of marching tetrahedra. """ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingTetrahedra, ::Type{VertType}=SVector{3,Float32}, ::Type{FaceType}=SVector{3, Int}; origin=VertType(-1,-1,-1), widths=VertType(2,2,2)) where {T, VertType, FaceType} - vts = Dict{Int, Int}() - fcs = FaceType[] + + vts = Dict{Int, Int}() + fcs = FaceType[] vtsAry = VertType[] sizehint!(vts, div(length(sdf),8)) sizehint!(vtsAry, div(length(sdf),8)) @@ -143,19 +143,23 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingTetrahedra, ::Type # process each voxel scale = widths ./ VertType(size(sdf) .- 1) nx::Int, ny::Int, nz::Int = size(sdf) + @inbounds for k = 1:nz-1, j = 1:ny-1, i= 1:nx-1 - vals = (sdf[i, j, k], - sdf[i, j+1, k], - sdf[i+1, j+1, k], - sdf[i+1, j, k], - sdf[i, j, k+1], - sdf[i, j+1, k+1], - sdf[i+1, j+1, k+1], - sdf[i+1, j, k+1]) - cubeindex = _get_cubeindex(vals,method.iso) - if cubeindex != 0x00 && cubeindex != 0xff - procVox(vals, method.iso, i, j, k, nx, ny, scale, origin, vts, vtsAry, fcs, method.eps, cubeindex) - end + + vals = (sdf[i, j, k ], + sdf[i, j+1,k ], + sdf[i+1,j+1,k ], + sdf[i+1,j, k ], + sdf[i, j, k+1], + sdf[i, j+1,k+1], + sdf[i+1,j+1,k+1], + sdf[i+1,j, k+1]) + + cubeindex = _get_cubeindex(vals, method.iso) + + _no_triangles(cubeindex) && continue + + procVox(vals, method.iso, i, j, k, nx, ny, scale, origin, vts, vtsAry, fcs, method.eps, cubeindex) end vtsAry,fcs From de4f2dfd5f9e27f1b7339fa024b5af8447916486 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Sat, 21 Dec 2019 22:56:27 -0500 Subject: [PATCH 17/19] [mt] support insidepositive sign conventions --- src/algorithmtypes.jl | 2 +- src/marching_tetrahedra.jl | 2 +- test/runtests.jl | 18 ++++++++++++++++++ 3 files changed, 20 insertions(+), 2 deletions(-) diff --git a/src/algorithmtypes.jl b/src/algorithmtypes.jl index d65174d..93680d1 100644 --- a/src/algorithmtypes.jl +++ b/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 diff --git a/src/marching_tetrahedra.jl b/src/marching_tetrahedra.jl index 907f76f..1c4722c 100644 --- a/src/marching_tetrahedra.jl +++ b/src/marching_tetrahedra.jl @@ -155,7 +155,7 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingTetrahedra, ::Type sdf[i+1,j+1,k+1], sdf[i+1,j, k+1]) - cubeindex = _get_cubeindex(vals, method.iso) + cubeindex = method.insidepositive ? _get_cubeindex_pos(vals, method.iso) : _get_cubeindex(vals, method.iso) _no_triangles(cubeindex) && continue diff --git a/test/runtests.jl b/test/runtests.jl index 30b69cd..8b7fa15 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -150,6 +150,24 @@ using LinearAlgebra: dot, norm @test length(faces(m)) == length(faces(mf)) end + @testset "sign flips" begin + algos = (MarchingCubes, ) #MarchingTetrahedra + + for algo_type in algos + algo = algo_type() + algo_pos = algo_type(insidepositive=true) + mf = SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)),algo, samples=(21,21,21)) do v + sqrt(sum(dot(v,v))) - 1 # sphere + end + + mf_pos = SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)),algo_pos, samples=(21,21,21)) do v + -sqrt(sum(dot(v,v))) + 1 # sphere positive inside + end + @test mf.vertices == mf_pos.vertices + @test mf.faces == mf_pos.faces + end + end + @testset "respect origin" begin # verify that when we construct a mesh, that mesh: # a) respects the origin of the SDF From 81fd819e1fde62c79095f06ba343b61be728dcef Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Sat, 21 Dec 2019 23:25:25 -0500 Subject: [PATCH 18/19] [mt] support direct meshing from function --- bench/benchmark.jl | 2 +- src/marching_tetrahedra.jl | 58 +++++++++++++++++++++++++++++++++++++- test/runtests.jl | 2 +- 3 files changed, 59 insertions(+), 3 deletions(-) diff --git a/bench/benchmark.jl b/bench/benchmark.jl index 81c1689..0ab9143 100644 --- a/bench/benchmark.jl +++ b/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 diff --git a/src/marching_tetrahedra.jl b/src/marching_tetrahedra.jl index 1c4722c..a61b1e4 100644 --- a/src/marching_tetrahedra.jl +++ b/src/marching_tetrahedra.jl @@ -144,7 +144,7 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingTetrahedra, ::Type scale = widths ./ VertType(size(sdf) .- 1) nx::Int, ny::Int, nz::Int = size(sdf) - @inbounds for k = 1:nz-1, j = 1:ny-1, i= 1:nx-1 + @inbounds for k = 1:nz-1, j = 1:ny-1, i = 1:nx-1 vals = (sdf[i, j, k ], sdf[i, j+1,k ], @@ -164,3 +164,59 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingTetrahedra, ::Type vtsAry,fcs end + +function isosurface(f::Function, method::MarchingTetrahedra, + ::Type{VertType}=SVector{3,Float32}, ::Type{FaceType}=SVector{3, Int}; + origin=VertType(-1,-1,-1), widths=VertType(2,2,2), + samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {T, VertType, FaceType} + + vts = Dict{Int, Int}() + fcs = FaceType[] + vtsAry = VertType[] + sizehint!(vts, div(prod(samples),8)) + sizehint!(vtsAry, div(prod(samples),8)) + sizehint!(fcs, div(prod(samples),4)) + # process each voxel + scale = widths ./ (VertType(samples...) .-1) + nx::Int, ny::Int, nz::Int = samples + zv = zero(eltype(VertType)) + vals = (zv,zv,zv,zv,zv,zv,zv,zv) + + @inbounds for k = 1:nz-1, j = 1:ny-1, i = 1:nx-1 + points = (VertType(i-1,j-1,k-1).* scale + origin, + VertType(i-1,j ,k-1).* scale + origin, + VertType(i ,j ,k-1).* scale + origin, + VertType(i ,j-1,k-1).* scale + origin, + VertType(i-1,j-1,k ).* scale + origin, + VertType(i-1,j ,k ).* scale + origin, + VertType(i ,j ,k ).* scale + origin, + VertType(i ,j-1,k ).* scale + origin) + if i == 0 + 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 + vals = (vals[4], + vals[3], + f(points[3]), + f(points[4]), + vals[8], + vals[7], + f(points[7]), + f(points[8])) + end + + cubeindex = method.insidepositive ? _get_cubeindex_pos(vals, method.iso) : _get_cubeindex(vals, method.iso) + + _no_triangles(cubeindex) && continue + + procVox(vals, method.iso, i, j, k, nx, ny, scale, origin, vts, vtsAry, fcs, method.eps, cubeindex) + end + + vtsAry,fcs +end diff --git a/test/runtests.jl b/test/runtests.jl index 8b7fa15..7b36048 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -151,7 +151,7 @@ using LinearAlgebra: dot, norm end @testset "sign flips" begin - algos = (MarchingCubes, ) #MarchingTetrahedra + algos = (MarchingCubes, MarchingTetrahedra) for algo_type in algos algo = algo_type() From cfc14f2d906373c2ee01e9983f2e74c4eddbb712 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Sun, 22 Dec 2019 11:16:12 -0500 Subject: [PATCH 19/19] [sn] support insidepositive --- docs/src/api.md | 3 ++- src/algorithmtypes.jl | 2 +- src/surface_nets.jl | 4 ++-- test/runtests.jl | 2 +- 4 files changed, 6 insertions(+), 5 deletions(-) diff --git a/docs/src/api.md b/docs/src/api.md index 13b5fa4..1ee11b7 100644 --- a/docs/src/api.md +++ b/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. diff --git a/src/algorithmtypes.jl b/src/algorithmtypes.jl index 93680d1..d21a58b 100644 --- a/src/algorithmtypes.jl +++ b/src/algorithmtypes.jl @@ -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 diff --git a/src/surface_nets.jl b/src/surface_nets.jl index 7a422b0..eb1b5fe 100644 --- a/src/surface_nets.jl +++ b/src/surface_nets.jl @@ -59,7 +59,7 @@ function isosurface(sdf::AbstractArray{T, 3}, method::NaiveSurfaceNets, ::Type{V sdf[xi+1,yi+2,zi+2], sdf[xi+2,yi+2,zi+2]) - mask = _get_cubeindex(grid, method.iso) + mask = method.insidepositive ? _get_cubeindex_pos(grid, method.iso) : _get_cubeindex(grid, method.iso) # Check for early termination if cell does not intersect boundary if _no_triangles(mask) @@ -168,7 +168,7 @@ function isosurface(f::Function, method::NaiveSurfaceNets, end # Also calculate 8-bit mask, like in marching cubes, so we can speed up sign checks later - mask = _get_cubeindex(grid, method.iso) + mask = method.insidepositive ? _get_cubeindex_pos(grid, method.iso) : _get_cubeindex(grid, method.iso) # Check for early termination if cell does not intersect boundary if _no_triangles(mask) diff --git a/test/runtests.jl b/test/runtests.jl index 7b36048..63af2f3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -151,7 +151,7 @@ using LinearAlgebra: dot, norm end @testset "sign flips" begin - algos = (MarchingCubes, MarchingTetrahedra) + algos = (MarchingCubes, MarchingTetrahedra, NaiveSurfaceNets) for algo_type in algos algo = algo_type()