From 7322fc22d79cf450fddc14f30715110f020ebc37 Mon Sep 17 00:00:00 2001
From: Steve Kelly <kd2cca@gmail.com>
Date: Fri, 14 Jun 2019 18:55:43 -0400
Subject: [PATCH 1/5] [MC] First pass at allocation-free meshing for MC

---
 src/marching_cubes.jl | 142 ++++++++++++++++++++++++++++++++++++++++++
 test/runtests.jl      |   9 +++
 2 files changed, 151 insertions(+)

diff --git a/src/marching_cubes.jl b/src/marching_cubes.jl
index e0df81a..fb88817 100644
--- a/src/marching_cubes.jl
+++ b/src/marching_cubes.jl
@@ -392,6 +392,144 @@ function marching_cubes(sdf::SignedDistanceField{3,ST,FT},
     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
+        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
+
+        # 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
+    end
+    MT(vts,fcs)
+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 +553,7 @@ 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
diff --git a/test/runtests.jl b/test/runtests.jl
index d47cb9e..f410db0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -76,6 +76,15 @@ using LinearAlgebra: dot, norm
         @test m == m2
 
     end
+
+    @testset "marching cubes function" begin
+        m = marching_cubes(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)),(21,21,21)) do v
+            sqrt(sum(dot(v,v))) - 1 # sphere
+        end
+        @test length(vertices(m)) == 10968
+        @test length(faces(m)) == 3656
+    end
+
     @testset "respect origin" begin
         # verify that when we construct a mesh, that mesh:
         #   a) respects the origin of the SDF

From 39c7d23be0f1c9154f94541dec2b8d436bea3bdd Mon Sep 17 00:00:00 2001
From: Steve Kelly <kd2cca@gmail.com>
Date: Sat, 15 Jun 2019 16:43:44 -0400
Subject: [PATCH 2/5] [SN] non-sdf SN implementation

---
 src/surface_nets.jl | 184 ++++++++++++++++++++++++++++++++++++++++++++
 test/runtests.jl    |  53 ++++++++-----
 2 files changed, 216 insertions(+), 21 deletions(-)

diff --git a/src/surface_nets.jl b/src/surface_nets.jl
index afdfff3..549da1b 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]<dims[3]-1
+
+        # m is the pointer into the buffer we are going to use.
+        # This is slightly obtuse because javascript does not have good support for packed data structures, so we must use typed arrays :(
+        # The contents of the buffer will be the indices of the vertices on the previous x/y slice of the volume
+        m = 1 + (dims[1]+1) * (1 + buf_no * (dims[2]+1))
+
+        x[2]=0
+        @inbounds while x[2]<dims[2]-1
+
+            x[1]=0
+            @inbounds while x[1] < dims[1]-1
+
+                # Read in 8 field values around this vertex and store them in an array
+                points = (Point{3,Float64}(x[1],x[2],x[3]).* scale + origin,
+                          Point{3,Float64}(x[1]+1,x[2],x[3]).* scale + origin,
+                          Point{3,Float64}(x[1],x[2]+1,x[3]).* scale + origin,
+                          Point{3,Float64}(x[1]+1,x[2]+1,x[3]).* scale + origin,
+                          Point{3,Float64}(x[1],x[2],x[3]+1).* scale + origin,
+                          Point{3,Float64}(x[1]+1,x[2],x[3]+1).* scale + origin,
+                          Point{3,Float64}(x[1],x[2]+1,x[3]+1).* scale + origin,
+                          Point{3,Float64}(x[1]+1,x[2]+1,x[3]+1).* scale + origin)
+
+                @inbounds grid = map(f, points)
+
+                # Also calculate 8-bit mask, like in marching cubes, so we can speed up sign checks later
+                mask = 0x00
+                signbit(grid[1]) && (mask |= 0x01)
+                signbit(grid[2]) && (mask |= 0x02)
+                signbit(grid[3]) && (mask |= 0x04)
+                signbit(grid[4]) && (mask |= 0x08)
+                signbit(grid[5]) && (mask |= 0x10)
+                signbit(grid[6]) && (mask |= 0x20)
+                signbit(grid[7]) && (mask |= 0x40)
+                signbit(grid[8]) && (mask |= 0x80)
+
+                # Check for early termination if cell does not intersect boundary
+                if mask == 0x00 || mask == 0xff
+                    x[1] += 1
+                    n += 1
+                    m += 1
+                    continue
+                end
+
+                #Sum up edge intersections
+                edge_mask = sn_edge_table[mask+1]
+                v[1] = 0.0
+                v[2] = 0.0
+                v[3] = 0.0
+                e_count = 0
+
+                #For every edge of the cube...
+                @inbounds for i=0:11
+
+                    #Use edge mask to check if it is crossed
+                    if (edge_mask & (1<<i)) == 0
+                        continue
+                    end
+
+                    #If it did, increment number of edge crossings
+                    e_count += 1
+
+                    #Now find the point of intersection
+                    e0 = cube_edges[(i<<1)+1]       #Unpack vertices
+                    e1 = cube_edges[(i<<1)+2]
+                    g0 = grid[e0+1]                 #Unpack grid values
+                    g1 = grid[e1+1]
+                    t  = g0 - g1                 #Compute point of intersection
+                    if abs(t) > 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<<i)) == 0
+                        continue
+                    end
+
+                    # i = axes we are point along.  iu, iv = orthogonal axes
+                    iu = (i+1)%3
+                    iv = (i+2)%3
+
+                    #If we are on a boundary, skip it
+                    if (x[iu+1] == 0 || x[iv+1] == 0)
+                        continue
+                    end
+
+                    #Otherwise, look up adjacent edges in buffer
+                    du = R[iu+1]
+                    dv = R[iv+1]
+
+                    #Remember to flip orientation depending on the sign of the corner.
+                    if (mask & 0x01) != 0x00
+                        push!(faces,Face{4,Int}(buffer[m+1]+1, buffer[m-du+1]+1, buffer[m-du-dv+1]+1, buffer[m-dv+1]+1));
+                    else
+                        push!(faces,Face{4,Int}(buffer[m+1]+1, buffer[m-dv+1]+1, buffer[m-du-dv+1]+1, buffer[m-du+1]+1));
+                    end
+                end
+                x[1] += 1
+                n += 1
+                m += 1
+            end
+            x[2] += 1
+            n += 1
+            m += 2
+        end
+        x[3] += 1
+        n+=dims[1]
+        buf_no = xor(buf_no,1)
+        R[3]=-R[3]
+    end
+    #All done!  Return the result
+
+    vertices, faces # faces are quads, indexed to vertices
+end
+
+
 struct NaiveSurfaceNets{T} <: AbstractMeshingAlgorithm
     iso::T
     eps::T
@@ -235,3 +404,18 @@ function (::Type{MT})(sdf::SignedDistanceField, method::NaiveSurfaceNets) where
                             orig)
     MT(vts, fcs)::MT
 end
+
+function (::Type{MT})(f::Function, bounds::HyperRectangle, size::NTuple{3,Int}, method::NaiveSurfaceNets) where {MT <: AbstractMesh}
+    orig = origin(bounds)
+    w = widths(bounds)
+    scale = w ./ Point(size .- 1)  # subtract 1 because an SDF with N points per side has N-1 cells
+
+    # TODO ISO val
+
+    vts, fcs = surface_nets(f,
+                            size,
+                            method.eps,
+                            scale,
+                            orig)
+    MT(vts, fcs)::MT
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index f410db0..1b0e95b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,18 +9,32 @@ using LinearAlgebra: dot, norm
 
 @testset "meshing" begin
     @testset "surface nets" begin
-          sdf_sphere = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v
-              sqrt(sum(dot(v,v))) - 1 # sphere
-          end
-          sdf_torus = SignedDistanceField(HyperRectangle(Vec(-2,-2,-2.),Vec(4,4,4.)), 0.05) do v
-              (sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25
-          end
-          sphere = HomogenousMesh(sdf_sphere, NaiveSurfaceNets())
-          torus = HomogenousMesh(sdf_torus, NaiveSurfaceNets())
-          @test length(vertices(sphere)) == 1832
-          @test length(vertices(torus)) == 5532
-          @test length(faces(sphere)) == 1830
-          @test length(faces(torus)) == 5532
+        sdf_sphere = SignedDistanceField(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.))) do v
+            sqrt(sum(dot(v,v))) - 1 # sphere
+        end
+        sdf_torus = SignedDistanceField(HyperRectangle(Vec(-2,-2,-2.),Vec(4,4,4.)), 0.05) do v
+            (sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25
+        end
+
+        snf_torus = HomogenousMesh(HyperRectangle(Vec(-2,-2,-2.),Vec(4,4,4.)), (81,81,81), NaiveSurfaceNets()) do v
+            (sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25
+        end
+
+        snf_sphere = HomogenousMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), (21,21,21), NaiveSurfaceNets()) do v
+            sqrt(sum(dot(v,v))) - 1 # sphere
+        end
+
+        sphere = HomogenousMesh(sdf_sphere, NaiveSurfaceNets())
+        torus = HomogenousMesh(sdf_torus, NaiveSurfaceNets())
+        @test length(vertices(sphere)) == 1832
+        @test length(vertices(torus)) == 5532
+        @test length(faces(sphere)) == 1830
+        @test length(faces(torus)) == 5532
+
+        @test length(vertices(sphere)) == length(vertices(snf_sphere))
+        @test length(vertices(torus)) == length(vertices(snf_torus))
+        @test length(faces(sphere)) == length(faces(snf_sphere))
+        @test length(faces(torus)) == length(faces(snf_torus))
     end
 
 
@@ -69,20 +83,17 @@ using LinearAlgebra: dot, norm
             sqrt(sum(dot(v,v))) - 1 # sphere
         end
 
+        mf = marching_cubes(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)),(21,21,21)) do v
+            sqrt(sum(dot(v,v))) - 1 # sphere
+        end
+
         m = marching_cubes(sdf,0)
         m2 = marching_cubes(sdf)
         @test length(vertices(m)) == 10968
         @test length(faces(m)) == 3656
         @test m == m2
-
-    end
-
-    @testset "marching cubes function" begin
-        m = marching_cubes(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)),(21,21,21)) do v
-            sqrt(sum(dot(v,v))) - 1 # sphere
-        end
-        @test length(vertices(m)) == 10968
-        @test length(faces(m)) == 3656
+        @test length(vertices(m)) == length(vertices(mf))
+        @test length(faces(m)) == length(faces(mf))
     end
 
     @testset "respect origin" begin

From 44741d7564793d7764604d92cfe4ce936af7986e Mon Sep 17 00:00:00 2001
From: steve <kd2cca@gmail.com>
Date: Sat, 22 Jun 2019 00:09:33 -0400
Subject: [PATCH 3/5] generalize vertex interpolation code to share between MC
 implementations

---
 src/marching_cubes.jl | 188 +++++++++++++++++-------------------------
 1 file changed, 76 insertions(+), 112 deletions(-)

diff --git a/src/marching_cubes.jl b/src/marching_cubes.jl
index fb88817..b044488 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,56 +313,32 @@ 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)
-        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
+        find_vertices_interp!(vertlist, points, iso_vals, cubeindex, iso, eps)
 
         # Create the triangle
         for i = 1:3:13
@@ -418,7 +379,6 @@ function marching_cubes(f::Function,
     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
@@ -468,54 +428,7 @@ function marching_cubes(f::Function,
         # 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
-        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
+        find_vertices_interp!(vertlist, points, iso_vals, cubeindex, iso, eps)
 
         # Create the triangle
         for i = 1:3:13
@@ -530,6 +443,57 @@ function marching_cubes(f::Function,
     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)

From 38ef674c5490addb84f07e3a2eb61e60757735d4 Mon Sep 17 00:00:00 2001
From: Steve Kelly <kd2cca@gmail.com>
Date: Sun, 23 Jun 2019 16:01:51 -0400
Subject: [PATCH 4/5] update README

---
 README.md | 22 ++++++++++++----------
 1 file changed, 12 insertions(+), 10 deletions(-)

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

From 70bc91ce875266baa7dcf46252ab5f5992264dcd Mon Sep 17 00:00:00 2001
From: Steve Kelly <kd2cca@gmail.com>
Date: Sun, 23 Jun 2019 16:02:33 -0400
Subject: [PATCH 5/5] add convience constructors with pre-specified sample size
 of 128

---
 src/marching_cubes.jl |  4 ++++
 src/surface_nets.jl   |  4 ++++
 test/runtests.jl      | 16 ++++++++++++++++
 3 files changed, 24 insertions(+)

diff --git a/src/marching_cubes.jl b/src/marching_cubes.jl
index b044488..54e1ecf 100644
--- a/src/marching_cubes.jl
+++ b/src/marching_cubes.jl
@@ -521,3 +521,7 @@ 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 549da1b..94dd4dc 100644
--- a/src/surface_nets.jl
+++ b/src/surface_nets.jl
@@ -419,3 +419,7 @@ function (::Type{MT})(f::Function, bounds::HyperRectangle, size::NTuple{3,Int},
                             orig)
     MT(vts, fcs)::MT
 end
+
+function (::Type{MT})(f::Function, bounds::HyperRectangle, method::NaiveSurfaceNets;size::NTuple{3,Int}=(128,128,128)) where {MT <: AbstractMesh}
+    MT(f,bounds,size,method)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 1b0e95b..e36573d 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -24,6 +24,14 @@ using LinearAlgebra: dot, norm
             sqrt(sum(dot(v,v))) - 1 # sphere
         end
 
+        # test convience constructors
+        HomogenousMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), NaiveSurfaceNets()) do v
+            sqrt(sum(dot(v,v))) - 1 # sphere
+        end
+        HomogenousMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), NaiveSurfaceNets(), size=(5,5,5)) do v
+            sqrt(sum(dot(v,v))) - 1 # sphere
+        end
+
         sphere = HomogenousMesh(sdf_sphere, NaiveSurfaceNets())
         torus = HomogenousMesh(sdf_torus, NaiveSurfaceNets())
         @test length(vertices(sphere)) == 1832
@@ -87,6 +95,14 @@ using LinearAlgebra: dot, norm
             sqrt(sum(dot(v,v))) - 1 # sphere
         end
 
+        # convience constructors
+        SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), MarchingCubes()) do v
+            sqrt(sum(dot(v,v))) - 1 # sphere
+        end
+        SimpleMesh(HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), MarchingCubes(), size=(5,6,7)) do v
+            sqrt(sum(dot(v,v))) - 1 # sphere
+        end
+
         m = marching_cubes(sdf,0)
         m2 = marching_cubes(sdf)
         @test length(vertices(m)) == 10968