diff --git a/etc/output_comparison.jl b/etc/output_comparison.jl
new file mode 100644
index 0000000..91b96fe
--- /dev/null
+++ b/etc/output_comparison.jl
@@ -0,0 +1,17 @@
+
+using Meshing
+using GeometryTypes
+using LinearAlgebra: dot, norm
+using FileIO
+
+f(v) = sqrt(sum(dot(v,v))) - 1
+sdf = SignedDistanceField(f,HyperRectangle(Vec(-1,-1,-1.), Vec(2,2,2.)))
+
+mc = HomogenousMesh(sdf, MarchingCubes())
+mt = HomogenousMesh(sdf, MarchingTetrahedra())
+ns = HomogenousMesh(sdf, NaiveSurfaceNets())
+
+# save the Sphere as a PLY file
+save("sphere_mc.ply",mc)
+save("sphere_mt.ply",mt)
+save("sphere_ns.ply",ns)
diff --git a/src/common.jl b/src/common.jl
index 1230e9a..039822b 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -3,7 +3,7 @@
     _get_cubeindex(iso_vals, iso)
 
 given `iso_vals` and iso, return an 8 bit value corresponding
-to each corner of a cube. In each bit position, 
+to each corner of a cube. In each bit position,
 0 indicates in the isosurface and 1 indicates outside the surface,
 where the sign convention indicates negative inside the surface
 """
@@ -23,7 +23,7 @@ end
     _get_cubeindex_pos(iso_vals, iso)
 
 given `iso_vals` and iso, return an 8 bit value corresponding
-to each corner of a cube. In each bit position, 
+to each corner of a cube. In each bit position,
 0 indicates in the isosurface and 1 indicates outside the surface,
 where the sign convention indicates positive inside the surface
 """
diff --git a/src/lut/sn.jl b/src/lut/sn.jl
index 041ad0b..55f8b4c 100644
--- a/src/lut/sn.jl
+++ b/src/lut/sn.jl
@@ -1,6 +1,5 @@
-const cube_edges = (1, 2, 1, 3, 1, 5, 2, 4, 2, 6, 3, 4,
-                    3, 7, 4, 8, 5, 6, 5, 7, 6, 8, 7, 8)
-
+const cube_edges = ( 0, 1, 0, 2, 0, 4, 1, 3, 1, 5, 2, 3,
+                     2, 6, 3, 7, 4, 5, 4, 6, 5, 7, 6, 7 )
 const sn_edge_table = (7, 25, 30, 98, 101, 123, 124, 168, 175, 177, 182, 202,
                        205, 211, 212, 772, 771, 797, 794, 870, 865, 895, 888,
                        940, 939, 949, 946, 974, 969, 983, 976, 1296, 1303, 1289,
diff --git a/src/surface_nets.jl b/src/surface_nets.jl
index 89a0a96..104414a 100644
--- a/src/surface_nets.jl
+++ b/src/surface_nets.jl
@@ -82,7 +82,7 @@ function isosurface(sdf::AbstractArray{T, 3}, method::NaiveSurfaceNets, ::Type{V
                 #Sum up edge intersections
                 edge_mask = sn_edge_table[mask]
 
-                _sn_add_verts!(xi, yi, zi, vertices, grid, edge_mask, buffer, m, scale, origin, method.eps, T, Val(true), VertType)
+                _sn_add_verts!(xi, yi, zi, vertices, grid, edge_mask, buffer, m, scale, origin, method.eps, Val(true), VertType)
 
                 #Now we need to add faces together, to do this we just loop over 3 basis components
                 x = (xi,yi,zi)
@@ -221,7 +221,7 @@ function isosurface(f::Function, method::NaiveSurfaceNets,
                 edge_mask = sn_edge_table[mask]
 
                 # add vertices
-                _sn_add_verts!(xi, yi, zi, vertices, grid, edge_mask, buffer, m, scale, origin, method.eps, eltype(VertType), Val(false), VertType)
+                _sn_add_verts!(xi, yi, zi, vertices, grid, edge_mask, buffer, m, scale, origin, method.eps, Val(false), VertType)
 
                 #Now we need to add faces together, to do this we just loop over 3 basis components
                 x = (xi,yi,zi)
@@ -269,8 +269,9 @@ function isosurface(f::Function, method::NaiveSurfaceNets,
     vertices, faces # faces are quads, indexed to vertices
 end
 
-@inline function _sn_add_verts!(xi, yi, zi, vertices, grid, edge_mask, buffer, m, scale, origin, eps, T, translate_pt, ::Type{VertType}) where {VertType}
+@inline function _sn_add_verts!(xi, yi, zi, vertices, grid, edge_mask, buffer, m, scale, origin, eps, translate_pt, ::Type{VertType}) where {VertType}
     v = zero(VertType)
+    T = eltype(VertType)
     e_count = 0
 
     #For every edge of the cube...
@@ -287,8 +288,8 @@ end
         #Now find the point of intersection
         e0 = cube_edges[(i<<1)+1]       #Unpack vertices
         e1 = cube_edges[(i<<1)+2]
-        g0 = grid[e0]                 #Unpack grid values
-        g1 = grid[e1]
+        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
@@ -303,25 +304,25 @@ end
         a = e0 & 1
         b = e1 & 1
         (a != 0) && (xj += one(T))
-        (a != b) && (xj += (a != 0 ? - t : t))
+        (a != b) && (xj += (a != 0 ? -t : t))
         a = e0 & 2
         b = e1 & 2
         (a != 0) && (yj += one(T))
-        (a != b) && (yj += (a != 0 ? - t : t))
+        (a != b) && (yj += (a != 0 ? -t : t))
         a = e0 & 4
         b = e1 & 4
-        (a != 0) && (zj += 1.0)
-        (a != b) && (zj += (a != 0 ? - t : t))
+        (a != 0) && (zj += one(T))
+        (a != b) && (zj += (a != 0 ? -t : t))
         v += VertType(xj,yj,zj)
 
     end # edge check
 
     #Now we just average the edge intersections and add them to coordinate
     s = 1.0 / e_count
-    if typeof(translate_pt) == Val{true}
-        @inbounds v = (VertType(xi,yi,zi)  + s .* v) .* scale + origin
+    if typeof(translate_pt) === Val{true}
+        v = (VertType(xi,yi,zi)  .+ s .* v) .* scale + origin
     else
-        @inbounds v = (VertType(xi,yi,zi) + s .* v)# * scale[i] + origin[i]
+        v = (VertType(xi,yi,zi) .+ s .* v)# * scale[i] + origin[i]
     end
 
     #Add vertex to buffer, store pointer to vertex index in buffer