Merge pull request #25 from JuliaGeometry/sjk/perf

[WIP] non-SDF Meshing

Author: sjkelly

README.md

@@ -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).
 

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