diff --git a/README.md b/README.md index 1c6b6f3..020e439 100644 --- a/README.md +++ b/README.md @@ -3,18 +3,19 @@ [![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 @@ -22,17 +23,18 @@ 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). diff --git a/src/marching_cubes.jl b/src/marching_cubes.jl index e0df81a..54e1ecf 100644 --- a/src/marching_cubes.jl +++ b/src/marching_cubes.jl @@ -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, @@ -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 @@ -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) @@ -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 \ No newline at end of file diff --git a/src/surface_nets.jl b/src/surface_nets.jl index afdfff3..94dd4dc 100644 --- a/src/surface_nets.jl +++ b/src/surface_nets.jl @@ -205,6 +205,175 @@ function surface_nets(data::Vector{T}, dims,eps,scale,origin) where {T} vertices, faces # faces are quads, indexed to vertices end +""" +Generate a mesh using naive surface nets. +This takes the center of mass of the voxel as the vertex for each cube. +""" +function surface_nets(f::Function, dims::NTuple{3,Int},eps,scale,origin) + + # TODO + T = Float64 + + vertices = Point{3,T}[] + faces = Face{4,Int}[] + + sizehint!(vertices,ceil(Int,maximum(dims)^2/2)) + sizehint!(faces,ceil(Int,maximum(dims)^2/2)) + + n = 0 + x = [0,0,0] + R = Array{Int}([1, (dims[1]+1), (dims[1]+1)*(dims[2]+1)]) + buf_no = 1 + + buffer = fill(zero(Int),R[3]*2) + + v = Vector{T}([0.0,0.0,0.0]) + + #March over the voxel grid + x[3] = 0 + @inbounds while x[3] eps + t = g0 / t + else + continue + end + + #Interpolate vertices and add up intersections (this can be done without multiplying) + k = 1 + for j = 1:3 + a = e0 & k + b = e1 & k + (a != 0) && (v[j] += 1.0) + if a != b + v[j] += (a != 0 ? - t : t) + end + k<<=1 + end + end # edge check + + #Now we just average the edge intersections and add them to coordinate + s = 1.0 / e_count + for i=1:3 + @inbounds v[i] = (x[i] + s * v[i])# * scale[i] + origin[i] + end + + #Add vertex to buffer, store pointer to vertex index in buffer + buffer[m+1] = length(vertices) + push!(vertices, Point{3,T}(v[1],v[2],v[3])) + + #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<