Skip to content

Sjk/cleanup1 #51

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 5 commits into from
Oct 18, 2019
Merged

Sjk/cleanup1 #51

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
132c40e
[mc] deduplicate some code for point generation
sjkelly Oct 18, 2019
611c0c5
[mt] cleanup some function signatures
sjkelly Oct 18, 2019
4fa7bfd
[mt] lut: uint8 is good for the environment
sjkelly Oct 18, 2019
59a981d
[mc] use tuples in voxel edge interpolation rather than array allocat…
sjkelly Oct 18, 2019
bd129a0
[mc] aesthetics
sjkelly Oct 18, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion docs/src/internals.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
```


Expand Down
110 changes: 55 additions & 55 deletions src/lut/mt.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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))
113 changes: 41 additions & 72 deletions src/marching_cubes.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand All @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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

"""
Expand All @@ -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
34 changes: 17 additions & 17 deletions src/marching_tetrahedra.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

"""
Expand All @@ -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
Expand All @@ -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.
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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

Expand Down Expand Up @@ -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

Expand Down