diff --git a/docs/src/internals.md b/docs/src/internals.md index 8ba8dcf..fa4e621 100644 --- a/docs/src/internals.md +++ b/docs/src/internals.md @@ -18,8 +18,9 @@ Meshing._determine_types ```@docs Meshing._mc_create_triangles! Meshing._mc_unique_triangles! -Meshing.find_vertices_interp! +Meshing.find_vertices_interp Meshing.vertex_interp +Meshing.mc_vert_points ``` diff --git a/src/lut/mt.jl b/src/lut/mt.jl index 48cb7aa..f787da7 100644 --- a/src/lut/mt.jl +++ b/src/lut/mt.jl @@ -33,72 +33,72 @@ Voxel corner and edge indexing conventions PT(1, 0, 1)) end -const voxCrnrPosInt = voxCrnrPos(SVector{3,Int}) +const voxCrnrPosInt = voxCrnrPos(SVector{3,UInt8}) # the voxel IDs at either end of the tetrahedra edges, by edge ID -const voxEdgeCrnrs = ((1, 2), - (2, 3), - (4, 3), - (1, 4), - (5, 6), - (6, 7), - (8, 7), - (5, 8), - (1, 5), - (2, 6), - (3, 7), - (4, 8), - (1, 3), - (1, 8), - (1, 6), - (5, 7), - (2, 7), - (4, 7), - (1, 7)) +const voxEdgeCrnrs = ((0x01, 0x02), + (0x02, 0x03), + (0x04, 0x03), + (0x01, 0x04), + (0x05, 0x06), + (0x06, 0x07), + (0x08, 0x07), + (0x05, 0x08), + (0x01, 0x05), + (0x02, 0x06), + (0x03, 0x07), + (0x04, 0x08), + (0x01, 0x03), + (0x01, 0x08), + (0x01, 0x06), + (0x05, 0x07), + (0x02, 0x07), + (0x04, 0x07), + (0x01, 0x07)) # direction codes: # 0 => +x, 1 => +y, 2 => +z, # 3 => +xy, 4 => +xz, 5 => +yz, 6 => +xyz -const voxEdgeDir = (1,0,1,0,1,0,1,0,2,2,2,2,3,4,5,3,4,5,6) +const voxEdgeDir = (0x01, 0x00, 0x01, 0x00, 0x01, 0x00, 0x01, 0x00, 0x02, 0x02, 0x02, 0x02, 0x03, 0x04, 0x05, 0x03, 0x04, 0x05, 0x06) # For a pair of corner IDs, the edge ID joining them # 0 denotes a pair with no edge -const voxEdgeIx = ((0,1,13,4,9,15,19,14), - (1,0,2,0,0,10,17,0), - (13,2,0,3,0,0,11,0), - (4,0,3,0,0,0,18,12), - (9,0,0,0,0,5,16,8), - (15,10,0,0,5,0,6,0), - (19,17,11,18,16,6,0,7), - (14,0,0,12,8,0,7,0)) +const voxEdgeIx = ((0x00, 0x01, 0x0d, 0x04, 0x09, 0x0f, 0x13, 0x0e), + (0x01, 0x00, 0x02, 0x00, 0x00, 0x0a, 0x11, 0x00), + (0x0d, 0x02, 0x00, 0x03, 0x00, 0x00, 0x0b, 0x00), + (0x04, 0x00, 0x03, 0x00, 0x00, 0x00, 0x12, 0x0c), + (0x09, 0x00, 0x00, 0x00, 0x00, 0x05, 0x10, 0x08), + (0x0f, 0x0a, 0x00, 0x00, 0x05, 0x00, 0x06, 0x00), + (0x13, 0x11, 0x0b, 0x12, 0x10, 0x06, 0x00, 0x07), + (0x0e, 0x00, 0x00, 0x0c, 0x08, 0x00, 0x07, 0x00)) # voxel corners that comprise each of the six tetrahedra -const subTets = ((1,3,2,7), - (1,8,4,7), - (1,4,3,7), - (1,2,6,7), - (1,5,8,7), - (1,6,5,7)) +const subTets = ((0x01, 0x03, 0x02, 0x07), + (0x01, 0x08, 0x04, 0x07), + (0x01, 0x04, 0x03, 0x07), + (0x01, 0x02, 0x06, 0x07), + (0x01, 0x05, 0x08, 0x07), + (0x01, 0x06, 0x05, 0x07)) # tetrahedron corners for each edge (indices 1-4) -const tetEdgeCrnrs = ((1,2), - (2,3), - (1,3), - (1,4), - (2,4), - (3,4)) +const tetEdgeCrnrs = ((0x01, 0x02), + (0x02, 0x03), + (0x01, 0x03), + (0x01, 0x04), + (0x02, 0x04), + (0x03, 0x04)) # triangle cases for a given tetrahedron edge code -const tetTri = ((1,3,4,0,0,0), - (1,5,2,0,0,0), - (3,5,2,3,4,5), - (2,6,3,0,0,0), - (1,6,4,1,2,6), - (1,5,6,1,6,3), - (4,5,6,0,0,0), - (4,6,5,0,0,0), - (1,6,5,1,3,6), - (1,4,6,1,6,2), - (2,3,6,0,0,0), - (3,2,5,3,5,4), - (1,2,5,0,0,0), - (1,4,3,0,0,0)) +const tetTri = ((0x01, 0x03, 0x04, 0x00, 0x00, 0x00), + (0x01, 0x05, 0x02, 0x00, 0x00, 0x00), + (0x03, 0x05, 0x02, 0x03, 0x04, 0x05), + (0x02, 0x06, 0x03, 0x00, 0x00, 0x00), + (0x01, 0x06, 0x04, 0x01, 0x02, 0x06), + (0x01, 0x05, 0x06, 0x01, 0x06, 0x03), + (0x04, 0x05, 0x06, 0x00, 0x00, 0x00), + (0x04, 0x06, 0x05, 0x00, 0x00, 0x00), + (0x01, 0x06, 0x05, 0x01, 0x03, 0x06), + (0x01, 0x04, 0x06, 0x01, 0x06, 0x02), + (0x02, 0x03, 0x06, 0x00, 0x00, 0x00), + (0x03, 0x02, 0x05, 0x03, 0x05, 0x04), + (0x01, 0x02, 0x05, 0x00, 0x00, 0x00), + (0x01, 0x04, 0x03, 0x00, 0x00, 0x00)) diff --git a/src/marching_cubes.jl b/src/marching_cubes.jl index fa553b7..654cb05 100644 --- a/src/marching_cubes.jl +++ b/src/marching_cubes.jl @@ -24,7 +24,6 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingCubes, ::Type{Vert method.reduceverts && sizehint!(vts, mt*mt*5) !method.reduceverts && sizehint!(vts, mt*mt*6) sizehint!(fcs, mt*mt*2) - vertlist = Vector{VertType}(undef, 12) @inbounds for zi = 1:nz-1, yi = 1:ny-1, xi = 1:nx-1 @@ -44,21 +43,14 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingCubes, ::Type{Vert # Cube is entirely in/out of the surface (cubeindex == 0x00 || cubeindex == 0xff) && continue - points = (VertType(xi-1,yi-1,zi-1) .* s .+ origin, - VertType(xi,yi-1,zi-1) .* s .+ origin, - VertType(xi,yi,zi-1) .* s .+ origin, - VertType(xi-1,yi,zi-1) .* s .+ origin, - VertType(xi-1,yi-1,zi) .* s .+ origin, - VertType(xi,yi-1,zi) .* s .+ origin, - VertType(xi,yi,zi) .* s .+ origin, - VertType(xi-1,yi,zi) .* s .+ origin) + points = mc_vert_points(xi,yi,zi,s,origin,VertType) # Find the vertices where the surface intersects the cube - find_vertices_interp!(vertlist, points, iso_vals, cubeindex, method.iso, method.eps) + vertlist = find_vertices_interp(points, iso_vals, cubeindex, method.iso, method.eps) # Create the triangle - method.reduceverts == true && _mc_unique_triangles!(vts, fcs, vertlist, cubeindex, FaceType) - method.reduceverts == false && _mc_create_triangles!(vts, fcs, vertlist, cubeindex, FaceType) + method.reduceverts && _mc_unique_triangles!(vts, fcs, vertlist, cubeindex, FaceType) + !method.reduceverts && _mc_create_triangles!(vts, fcs, vertlist, cubeindex, FaceType) end vts,fcs end @@ -80,18 +72,11 @@ 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) - vertlist = Vector{VertType}(undef, 12) iso_vals = Vector{eltype(VertType)}(undef,8) @inbounds for xi = 1:nx-1, yi = 1:ny-1, zi = 1:nz-1 - points = (VertType(xi-1,yi-1,zi-1) .* s .+ origin, - VertType(xi,yi-1,zi-1) .* s .+ origin, - VertType(xi,yi,zi-1) .* s .+ origin, - VertType(xi-1,yi,zi-1) .* s .+ origin, - VertType(xi-1,yi-1,zi) .* s .+ origin, - VertType(xi,yi-1,zi) .* s .+ origin, - VertType(xi,yi,zi) .* s .+ origin, - VertType(xi-1,yi,zi) .* s .+ origin) + + points = mc_vert_points(xi,yi,zi,s,origin,VertType) if zi == 1 for i = 1:8 @@ -118,7 +103,7 @@ function isosurface(f::Function, method::MarchingCubes, ::Type{VertType}=SVector # Find the vertices where the surface intersects the cube # TODO this can use the underlying function to find the zero. # The underlying space is non-linear so there will be error otherwise - find_vertices_interp!(vertlist, points, iso_vals, cubeindex, method.iso, method.eps) + vertlist = find_vertices_interp(points, iso_vals, cubeindex, method.iso, method.eps) # Create the triangle method.reduceverts && _mc_unique_triangles!(vts, fcs, vertlist, cubeindex, FaceType) @@ -211,59 +196,27 @@ Each face may share a reference to a vertex with another face. end """ - find_vertices_interp!(vertlist, points, iso_vals, cubeindex, iso, eps) + find_vertices_interp(points, iso_vals, cubeindex, iso, eps) Find the vertices where the surface intersects the cube """ -@inline function find_vertices_interp!(vertlist, points, iso_vals, cubeindex, iso, eps) - if !iszero(edge_table[cubeindex] & 0x001) - vertlist[1] = - vertex_interp(iso,points[1],points[2],iso_vals[1],iso_vals[2], eps) - end - if !iszero(edge_table[cubeindex] & 0x002) - vertlist[2] = - vertex_interp(iso,points[2],points[3],iso_vals[2],iso_vals[3], eps) - end - if !iszero(edge_table[cubeindex] & 0x004) - vertlist[3] = - vertex_interp(iso,points[3],points[4],iso_vals[3],iso_vals[4], eps) - end - if !iszero(edge_table[cubeindex] & 0x008) - vertlist[4] = - vertex_interp(iso,points[4],points[1],iso_vals[4],iso_vals[1], eps) - end - if !iszero(edge_table[cubeindex] & 0x010) - vertlist[5] = - vertex_interp(iso,points[5],points[6],iso_vals[5],iso_vals[6], eps) - end - if !iszero(edge_table[cubeindex] & 0x020) - vertlist[6] = - vertex_interp(iso,points[6],points[7],iso_vals[6],iso_vals[7], eps) - end - if !iszero(edge_table[cubeindex] & 0x040) - vertlist[7] = - vertex_interp(iso,points[7],points[8],iso_vals[7],iso_vals[8], eps) - end - if !iszero(edge_table[cubeindex] & 0x080) - vertlist[8] = - vertex_interp(iso,points[8],points[5],iso_vals[8],iso_vals[5], eps) - end - if !iszero(edge_table[cubeindex] & 0x100) - vertlist[9] = - vertex_interp(iso,points[1],points[5],iso_vals[1],iso_vals[5], eps) - end - if !iszero(edge_table[cubeindex] & 0x200) - vertlist[10] = - vertex_interp(iso,points[2],points[6],iso_vals[2],iso_vals[6], eps) - end - if !iszero(edge_table[cubeindex] & 0x400) - vertlist[11] = - vertex_interp(iso,points[3],points[7],iso_vals[3],iso_vals[7], eps) - end - if !iszero(edge_table[cubeindex] & 0x800) - vertlist[12] = - vertex_interp(iso,points[4],points[8],iso_vals[4],iso_vals[8], eps) - end +@inline function find_vertices_interp(points, iso_vals, cubeindex, iso, eps) + VT = eltype(points) + zv = Ref{VT}()[] + # since we don't check any non-interpolated entries for zero elements (tri_table tells us the entries to check) + # we can use undefined references which is faster than zeroing out the memory + (!iszero(edge_table[cubeindex] & 0x001) ? vertex_interp(iso,points[1],points[2],iso_vals[1],iso_vals[2], eps) : zv, + !iszero(edge_table[cubeindex] & 0x002) ? vertex_interp(iso,points[2],points[3],iso_vals[2],iso_vals[3], eps) : zv, + !iszero(edge_table[cubeindex] & 0x004) ? vertex_interp(iso,points[3],points[4],iso_vals[3],iso_vals[4], eps) : zv, + !iszero(edge_table[cubeindex] & 0x008) ? vertex_interp(iso,points[4],points[1],iso_vals[4],iso_vals[1], eps) : zv, + !iszero(edge_table[cubeindex] & 0x010) ? vertex_interp(iso,points[5],points[6],iso_vals[5],iso_vals[6], eps) : zv, + !iszero(edge_table[cubeindex] & 0x020) ? vertex_interp(iso,points[6],points[7],iso_vals[6],iso_vals[7], eps) : zv, + !iszero(edge_table[cubeindex] & 0x040) ? vertex_interp(iso,points[7],points[8],iso_vals[7],iso_vals[8], eps) : zv, + !iszero(edge_table[cubeindex] & 0x080) ? vertex_interp(iso,points[8],points[5],iso_vals[8],iso_vals[5], eps) : zv, + !iszero(edge_table[cubeindex] & 0x100) ? vertex_interp(iso,points[1],points[5],iso_vals[1],iso_vals[5], eps) : zv, + !iszero(edge_table[cubeindex] & 0x200) ? vertex_interp(iso,points[2],points[6],iso_vals[2],iso_vals[6], eps) : zv, + !iszero(edge_table[cubeindex] & 0x400) ? vertex_interp(iso,points[3],points[7],iso_vals[3],iso_vals[7], eps) : zv, + !iszero(edge_table[cubeindex] & 0x800) ? vertex_interp(iso,points[4],points[8],iso_vals[4],iso_vals[8], eps) : zv) end """ @@ -282,3 +235,19 @@ function vertex_interp(iso, p1, p2, valp1, valp2, eps = 0.00001) return p end + +""" + mc_vert_points(xi,yi,zi, scale, origin, ::Type{VertType}) + +Returns a tuple of 8 points corresponding to each corner of a cube +""" +@inline function mc_vert_points(xi,yi,zi, scale, origin, ::Type{VertType}) where VertType + (VertType(xi-1,yi-1,zi-1) .* scale .+ origin, + VertType(xi ,yi-1,zi-1) .* scale .+ origin, + VertType(xi ,yi ,zi-1) .* scale .+ origin, + VertType(xi-1,yi ,zi-1) .* scale .+ origin, + VertType(xi-1,yi-1,zi ) .* scale .+ origin, + 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 diff --git a/src/marching_tetrahedra.jl b/src/marching_tetrahedra.jl index 7b53262..b1848b1 100644 --- a/src/marching_tetrahedra.jl +++ b/src/marching_tetrahedra.jl @@ -31,7 +31,7 @@ for vertex sharing to be implemented properly. """ @inline function vertId(e, x, y, z, nx, ny) dx, dy, dz = voxCrnrPosInt[voxEdgeCrnrs[e][1]] - voxEdgeDir[e]+7*(x-1+dx+nx*(y-1+dy+ny*(z-1+dz))) + voxEdgeDir[e]+7*(x-0x01+dx+nx*(y-0x01+dy+ny*(z-0x01+dz))) end """ @@ -49,7 +49,7 @@ corners (thereby preventing degeneracies). tgtVal = vals[ixs[2]] a = min(max((iso-srcVal)/(tgtVal-srcVal), eps), one(T)-eps) b = one(T)-a - c1= voxCrnrPos(VertType)[ixs[1]] + c1 = voxCrnrPos(VertType)[ixs[1]] c2 = voxCrnrPos(VertType)[ixs[2]] VertType(x,y,z) + c1 .* b + c2.* a @@ -60,7 +60,7 @@ end vals, iso::Real, vts::Dict, vtsAry::Vector, - eps::Real, ::Type{VertType}) where {VertType} + eps::Real) Gets the vertex ID, adding it to the vertex dictionary if not already present. @@ -69,8 +69,9 @@ present. vals, iso::Real, vts::Dict, vtsAry::Vector, - eps::Real, ::Type{VertType}) where {VertType} + eps::Real) + 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) @@ -95,17 +96,16 @@ end """ procVox(vals, iso::Real, x, y, z, nx, ny, vts::Dict, vtsAry::Vector, fcs::Vector, - eps::Real, - ::Type{VertType}, ::Type{FaceType}) where {VertType, FaceType} + eps::Real) 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, vts::Dict, vtsAry::Vector, fcs::Vector, - eps::Real, - ::Type{VertType}, ::Type{FaceType}) where {VertType, FaceType} - + eps::Real) + VertType = eltype(vtsAry) + FaceType = eltype(fcs) # check each sub-tetrahedron in the voxel @inbounds for i = 1:6 tIx = tetIx(i, vals, iso) @@ -115,16 +115,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, VertType), - getVertId(voxEdgeId(i, e[2]), x, y, z, nx, ny, vals, iso, vts, vtsAry, eps, VertType), - getVertId(voxEdgeId(i, e[3]), x, y, z, nx, ny, vals, iso, vts, vtsAry, eps, VertType))) + 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))) # bail if there are no more faces - e[4] == 0 && continue + iszero(e[4]) && continue push!(fcs, FaceType( - getVertId(voxEdgeId(i, e[4]), x, y, z, nx, ny, vals, iso, vts, vtsAry, eps, VertType), - getVertId(voxEdgeId(i, e[5]), x, y, z, nx, ny, vals, iso, vts, vtsAry, eps, VertType), - getVertId(voxEdgeId(i, e[6]), x, y, z, nx, ny, vals, iso, vts, vtsAry, eps, VertType))) + 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))) end end @@ -154,7 +154,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, VertType, FaceType) + procVox(vals, method.iso, i, j, k, nx, ny, vts, vtsAry, fcs, method.eps) end end