Skip to content

Commit

Permalink
Merge pull request #25 from JuliaGeometry/sjk/perf
Browse files Browse the repository at this point in the history 
[WIP] non-SDF Meshing
  • Loading branch information
@sjkelly
sjkelly authored Jun 23, 2019
2 parents 1c1a3ee + 70bc91c commit f1412a4
Show file tree
Hide file tree
Showing 4 changed files with 420 additions and 84 deletions.
22 changes: 12 additions & 10 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,36 +3,38 @@
[![Build Status](https://travis-ci.org/JuliaGeometry/Meshing.jl.svg)](https://travis-ci.org/JuliaGeometry/Meshing.jl)
[![codecov.io](http://codecov.io/github/JuliaGeometry/Meshing.jl/coverage.svg?branch=master)](http://codecov.io/github/JuliaGeometry/Meshing.jl?branch=master)

This package provides meshing algorithms for use on distance fields.
This package provides a comprehensive suite of meshing algorithms for use on distance fields.

Including:
Algorithms included:
* [Marching Tetrahedra](https://en.wikipedia.org/wiki/Marching_tetrahedra)
* [Marching Cubes](https://en.wikipedia.org/wiki/Marching_cubes)
* [Naive Surface Nets](https://0fps.net/2012/07/12/smooth-voxel-terrain-part-2/)

## Interface

This package is tightly integrated with [GeometryTypes.jl](https://github.com/JuliaGeometry/GeometryTypes.jl).
This package inherits the [mesh types](http://juliageometry.github.io/GeometryTypes.jl/latest/types.html#Meshes-1)
from [GeometryTypes.jl](https://github.com/JuliaGeometry/GeometryTypes.jl).

All algorithms operate on `SignedDistanceField` and output a concrete `AbstractMesh`. For example:
The algorithms operate on a `Function` or a `SignedDistanceField` and output a concrete `AbstractMesh`. For example:

```
using Meshing
using GeometryTypes
using LinearAlgebra: dot, norm
using FileIO
# generate an SDF of a sphere
sdf_sphere = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v
sqrt(sum(dot(v,v))) - 1 # sphere
# Mesh an equation of sphere in the Axis-Aligned Bounding box starting
# at -1,-1,-1 and widths of 2,2,2
m = GLNormalMesh(HyperRectangle(Vec(-1,-1,-1.), Vec(2,2,2.)), MarchingCubes()) do v
sqrt(sum(dot(v,v))) - 1
end
m = GLNormalMesh(sdf_sphere, MarchingCubes())
# save the Sphere as a PLY file
save("sphere.ply",m)
```

The general API is ``(::Type{MT})(sdf::SignedDistanceField, method::AbstractMeshingAlgorithm) where {MT <: AbstractMesh}``
The general API is: ```(::Type{MT})(sdf::Function, method::AbstractMeshingAlgorithm) where {MT <: AbstractMesh}``` or ```(::Type{MT})(sdf::SignedDistanceField, method::AbstractMeshingAlgorithm) where {MT <: AbstractMesh}```


For a full listing of concrete `AbstractMesh` types see [GeometryTypes.jl mesh documentation](http://juliageometry.github.io/GeometryTypes.jl/latest/types.html#Meshes-1).

Expand Down
232 changes: 171 additions & 61 deletions src/marching_cubes.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -305,21 +305,6 @@ function marching_cubes(sdf::SignedDistanceField{3,ST,FT},
vertlist = Vector{Point{3,Float64}}(undef, 12)
@inbounds for xi = 1:nx-1, yi = 1:ny-1, zi = 1:nz-1

#Determine the index into the edge table which
#tells us which vertices are inside of the surface
cubeindex = sdf[xi,yi,zi] < iso ? 1 : 0
sdf[xi+1,yi,zi] < iso && (cubeindex |= 2)
sdf[xi+1,yi+1,zi] < iso && (cubeindex |= 4)
sdf[xi,yi+1,zi] < iso && (cubeindex |= 8)
sdf[xi,yi,zi+1] < iso && (cubeindex |= 16)
sdf[xi+1,yi,zi+1] < iso && (cubeindex |= 32)
sdf[xi+1,yi+1,zi+1] < iso && (cubeindex |= 64)
sdf[xi,yi+1,zi+1] < iso && (cubeindex |= 128)
cubeindex += 1

# Cube is entirely in/out of the surface
edge_table[cubeindex] == 0 && continue

points = (Point{3,Float64}(xi-1,yi-1,zi-1) .* s .+ orig,
Point{3,Float64}(xi,yi-1,zi-1) .* s .+ orig,
Point{3,Float64}(xi,yi,zi-1) .* s .+ orig,
Expand All @@ -328,57 +313,123 @@ function marching_cubes(sdf::SignedDistanceField{3,ST,FT},
Point{3,Float64}(xi,yi-1,zi) .* s .+ orig,
Point{3,Float64}(xi,yi,zi) .* s .+ orig,
Point{3,Float64}(xi-1,yi,zi) .* s .+ orig)
iso_vals = (sdf[xi,yi,zi],
sdf[xi+1,yi,zi],
sdf[xi+1,yi+1,zi],
sdf[xi,yi+1,zi],
sdf[xi,yi,zi+1],
sdf[xi+1,yi,zi+1],
sdf[xi+1,yi+1,zi+1],
sdf[xi,yi+1,zi+1])

#Determine the index into the edge table which
#tells us which vertices are inside of the surface
cubeindex = iso_vals[1] < iso ? 1 : 0
iso_vals[2] < iso && (cubeindex |= 2)
iso_vals[3] < iso && (cubeindex |= 4)
iso_vals[4] < iso && (cubeindex |= 8)
iso_vals[5] < iso && (cubeindex |= 16)
iso_vals[6] < iso && (cubeindex |= 32)
iso_vals[7] < iso && (cubeindex |= 64)
iso_vals[8] < iso && (cubeindex |= 128)
cubeindex += 1

# Cube is entirely in/out of the surface
edge_table[cubeindex] == 0 && continue

# Find the vertices where the surface intersects the cube
if (edge_table[cubeindex] & 1 != 0)
vertlist[1] =
vertex_interp(iso,points[1],points[2],sdf[xi,yi,zi],sdf[xi+1,yi,zi], eps)
end
if (edge_table[cubeindex] & 2 != 0)
vertlist[2] =
vertex_interp(iso,points[2],points[3],sdf[xi+1,yi,zi],sdf[xi+1,yi+1,zi], eps)
end
if (edge_table[cubeindex] & 4 != 0)
vertlist[3] =
vertex_interp(iso,points[3],points[4],sdf[xi+1,yi+1,zi],sdf[xi,yi+1,zi], eps)
end
if (edge_table[cubeindex] & 8 != 0)
vertlist[4] =
vertex_interp(iso,points[4],points[1],sdf[xi,yi+1,zi],sdf[xi,yi,zi], eps)
end
if (edge_table[cubeindex] & 16 != 0)
vertlist[5] =
vertex_interp(iso,points[5],points[6],sdf[xi,yi,zi+1],sdf[xi+1,yi,zi+1], eps)
end
if (edge_table[cubeindex] & 32 != 0)
vertlist[6] =
vertex_interp(iso,points[6],points[7],sdf[xi+1,yi,zi+1],sdf[xi+1,yi+1,zi+1], eps)
end
if (edge_table[cubeindex] & 64 != 0)
vertlist[7] =
vertex_interp(iso,points[7],points[8],sdf[xi+1,yi+1,zi+1],sdf[xi,yi+1,zi+1], eps)
end
if (edge_table[cubeindex] & 128 != 0)
vertlist[8] =
vertex_interp(iso,points[8],points[5],sdf[xi,yi+1,zi+1],sdf[xi,yi,zi+1], eps)
end
if (edge_table[cubeindex] & 256 != 0)
vertlist[9] =
vertex_interp(iso,points[1],points[5],sdf[xi,yi,zi],sdf[xi,yi,zi+1], eps)
end
if (edge_table[cubeindex] & 512 != 0)
vertlist[10] =
vertex_interp(iso,points[2],points[6],sdf[xi+1,yi,zi],sdf[xi+1,yi,zi+1], eps)
end
if (edge_table[cubeindex] & 1024 != 0)
vertlist[11] =
vertex_interp(iso,points[3],points[7],sdf[xi+1,yi+1,zi],sdf[xi+1,yi+1,zi+1], eps)
find_vertices_interp!(vertlist, points, iso_vals, cubeindex, iso, eps)

# Create the triangle
for i = 1:3:13
tri_table[cubeindex][i] == -1 && break
push!(vts, vertlist[tri_table[cubeindex][i ]])
push!(vts, vertlist[tri_table[cubeindex][i+1]])
push!(vts, vertlist[tri_table[cubeindex][i+2]])
fct = length(vts)
push!(fcs, Face{3,Int}(fct, fct-1, fct-2))
end
if (edge_table[cubeindex] & 2048 != 0)
vertlist[12] =
vertex_interp(iso,points[4],points[8],sdf[xi,yi+1,zi],sdf[xi,yi+1,zi+1], eps)
end
MT(vts,fcs)
end


function marching_cubes(f::Function,
bounds::HyperRectangle,
samples::NTuple{3,Int}=(256,256,256),
iso=0.0,
MT::Type{M}=SimpleMesh{Point{3,Float64},Face{3,Int}},
eps=0.00001) where {ST,FT,M<:AbstractMesh}
nx, ny, nz = samples[1], samples[2], samples[3]
w = widths(bounds)
orig = origin(bounds)

# we subtract one from the length along each axis because
# an NxNxN SDF has N-1 cells on each axis
s = Point{3,Float64}(w[1]/(nx-1), w[2]/(ny-1), w[3]/(nz-1))

# arrays for vertices and faces
vts = Point{3,Float64}[]
fcs = Face{3,Int}[]
mt = max(nx,ny,nz)
sizehint!(vts, mt*mt*6)
sizehint!(fcs, mt*mt*2)
vertlist = Vector{Point{3,Float64}}(undef, 12)
iso_vals = Vector{Float64}(undef,8)
points = Vector{Point{3,Float64}}(undef,8)
@inbounds for xi = 1:nx-1, yi = 1:ny-1, zi = 1:nz-1

if zi == 1
points[1] = Point{3,Float64}(xi-1,yi-1,zi-1) .* s .+ orig
points[2] = Point{3,Float64}(xi,yi-1,zi-1) .* s .+ orig
points[3] = Point{3,Float64}(xi,yi,zi-1) .* s .+ orig
points[4] = Point{3,Float64}(xi-1,yi,zi-1) .* s .+ orig
points[5] = Point{3,Float64}(xi-1,yi-1,zi) .* s .+ orig
points[6] = Point{3,Float64}(xi,yi-1,zi) .* s .+ orig
points[7] = Point{3,Float64}(xi,yi,zi) .* s .+ orig
points[8] = Point{3,Float64}(xi-1,yi,zi) .* s .+ orig
for i = 1:8
iso_vals[i] = f(points[i])
end
else
points[1] = points[5]
points[2] = points[6]
points[3] = points[7]
points[4] = points[8]
points[5] = Point{3,Float64}(xi-1,yi-1,zi) .* s .+ orig
points[6] = Point{3,Float64}(xi,yi-1,zi) .* s .+ orig
points[7] = Point{3,Float64}(xi,yi,zi) .* s .+ orig
points[8] = Point{3,Float64}(xi-1,yi,zi) .* s .+ orig
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])
end

#Determine the index into the edge table which
#tells us which vertices are inside of the surface
cubeindex = iso_vals[1] < iso ? 1 : 0
iso_vals[2] < iso && (cubeindex |= 2)
iso_vals[3] < iso && (cubeindex |= 4)
iso_vals[4] < iso && (cubeindex |= 8)
iso_vals[5] < iso && (cubeindex |= 16)
iso_vals[6] < iso && (cubeindex |= 32)
iso_vals[7] < iso && (cubeindex |= 64)
iso_vals[8] < iso && (cubeindex |= 128)
cubeindex += 1

# Cube is entirely in/out of the surface
edge_table[cubeindex] == 0 && continue

# 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, iso, eps)

# Create the triangle
for i = 1:3:13
tri_table[cubeindex][i] == -1 && break
Expand All @@ -392,6 +443,57 @@ function marching_cubes(sdf::SignedDistanceField{3,ST,FT},
MT(vts,fcs)
end

@inline function find_vertices_interp!(vertlist, points, iso_vals, cubeindex, iso, eps)
if (edge_table[cubeindex] & 1 != 0)
vertlist[1] =
vertex_interp(iso,points[1],points[2],iso_vals[1],iso_vals[2], eps)
end
if (edge_table[cubeindex] & 2 != 0)
vertlist[2] =
vertex_interp(iso,points[2],points[3],iso_vals[2],iso_vals[3], eps)
end
if (edge_table[cubeindex] & 4 != 0)
vertlist[3] =
vertex_interp(iso,points[3],points[4],iso_vals[3],iso_vals[4], eps)
end
if (edge_table[cubeindex] & 8 != 0)
vertlist[4] =
vertex_interp(iso,points[4],points[1],iso_vals[4],iso_vals[1], eps)
end
if (edge_table[cubeindex] & 16 != 0)
vertlist[5] =
vertex_interp(iso,points[5],points[6],iso_vals[5],iso_vals[6], eps)
end
if (edge_table[cubeindex] & 32 != 0)
vertlist[6] =
vertex_interp(iso,points[6],points[7],iso_vals[6],iso_vals[7], eps)
end
if (edge_table[cubeindex] & 64 != 0)
vertlist[7] =
vertex_interp(iso,points[7],points[8],iso_vals[7],iso_vals[8], eps)
end
if (edge_table[cubeindex] & 128 != 0)
vertlist[8] =
vertex_interp(iso,points[8],points[5],iso_vals[8],iso_vals[5], eps)
end
if (edge_table[cubeindex] & 256 != 0)
vertlist[9] =
vertex_interp(iso,points[1],points[5],iso_vals[1],iso_vals[5], eps)
end
if (edge_table[cubeindex] & 512 != 0)
vertlist[10] =
vertex_interp(iso,points[2],points[6],iso_vals[2],iso_vals[6], eps)
end
if (edge_table[cubeindex] & 1024 != 0)
vertlist[11] =
vertex_interp(iso,points[3],points[7],iso_vals[3],iso_vals[7], eps)
end
if (edge_table[cubeindex] & 2048 != 0)
vertlist[12] =
vertex_interp(iso,points[4],points[8],iso_vals[4],iso_vals[8], eps)
end
end

# Linearly interpolate the position where an isosurface cuts
# an edge between two vertices, each with their own scalar value
function vertex_interp(iso, p1, p2, valp1, valp2, eps = 0.00001)
Expand All @@ -415,3 +517,11 @@ MarchingCubes(iso::T1=0.0, eps::T2=1e-3) where {T1, T2} = MarchingCubes{promote_
function (::Type{MT})(df::SignedDistanceField, method::MarchingCubes)::MT where {MT <: AbstractMesh}
marching_cubes(df, method.iso, MT, method.eps)
end

function (::Type{MT})(f::Function, h::HyperRectangle, size::NTuple{3,Int}, method::MarchingCubes)::MT where {MT <: AbstractMesh}
marching_cubes(f, h, size, method.iso, MT, method.eps)
end

function (::Type{MT})(f::Function, h::HyperRectangle, method::MarchingCubes; size::NTuple{3,Int}=(128,128,128))::MT where {MT <: AbstractMesh}
marching_cubes(f, h, size, method.iso, MT, method.eps)
end
Loading

0 comments on commit f1412a4

Please sign in to comment.