diff --git a/Project.toml b/Project.toml
index b44c442..6faa49a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -3,10 +3,8 @@ uuid = "e6723b4c-ebff-59f1-b4b7-d97aa5274f73"
 version = "0.7.0"
 
 [deps]
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
-StaticArrays = "1"
 julia = "1.9"
 
 [extras]
diff --git a/docs/src/internals.md b/docs/src/internals.md
index 572bbf2..4bfedca 100644
--- a/docs/src/internals.md
+++ b/docs/src/internals.md
@@ -8,8 +8,6 @@ This is some brief documentation on the basic internal functions.
 ```@docs
 Meshing._DEFAULT_SAMPLES
 Meshing._get_cubeindex
-Meshing._get_cubeindex_pos
-Meshing._determine_types
 ```
 
 ## Marching Cubes
diff --git a/src/Meshing.jl b/src/Meshing.jl
index eb3ff95..dbde984 100644
--- a/src/Meshing.jl
+++ b/src/Meshing.jl
@@ -1,7 +1,5 @@
 module Meshing
 
-using StaticArrays
-
 """
     _DEFAULT_SAMPLES = (24,24,24)
 
@@ -13,13 +11,11 @@ include("algorithmtypes.jl")
 include("common.jl")
 include("marching_tetrahedra.jl")
 include("marching_cubes.jl")
-include("surface_nets.jl")
 include("roots.jl")
 include("adaptive.jl")
 
 export isosurface,
        MarchingCubes,
-       MarchingTetrahedra,
-       NaiveSurfaceNets
+       MarchingTetrahedra
 
 end # module
diff --git a/src/adaptive.jl b/src/adaptive.jl
index 29b3afe..f2e6619 100644
--- a/src/adaptive.jl
+++ b/src/adaptive.jl
@@ -12,7 +12,7 @@ end
 
 
 function vertices(h::HyperRectangle, ::Type{SV}) where SV
-    z = zero(eltype(SV))
+    z = zero(eltype(h.origin))
     @inbounds begin
         o = SV(h.origin[1],h.origin[2],h.origin[3])
         w = SV(h.widths[1],h.widths[2],h.widths[3])
diff --git a/src/algorithmtypes.jl b/src/algorithmtypes.jl
index d1b81dc..02d0bf3 100644
--- a/src/algorithmtypes.jl
+++ b/src/algorithmtypes.jl
@@ -13,27 +13,8 @@ abstract type AbstractMeshingAlgorithm end
 abstract type AbstractAdaptiveMeshingAlgorithm end
 
 
-function (::Type{MeshAlgo})(;iso::T1=0.0, eps::T2=1e-3, reduceverts::Bool=true, insidepositive::Bool=false) where {T1, T2, MeshAlgo <: AbstractMeshingAlgorithm}
-    if isconcretetype(MeshAlgo)
-        return MeshAlgo(iso, eps, reduceverts, insidepositive)
-    else
-        return MeshAlgo{promote_type(T1, T2)}(iso, eps, reduceverts, insidepositive)
-    end
-end
-
-function (::Type{MeshAlgo})(iso) where MeshAlgo <: AbstractMeshingAlgorithm
-    MeshAlgo(iso=iso)
-end
-
-function (::Type{MeshAlgo})(iso,eps) where MeshAlgo <: AbstractMeshingAlgorithm
-    MeshAlgo(iso=iso,eps=eps)
-end
-
 """
-    MarchingCubes(iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
-    MarchingCubes(;iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
-    MarchingCubes(iso)
-    MarchingCubes(iso,eps)
+    MarchingCubes(;iso=0.0)
 
 Specifies the use of the Marching Cubes algorithm for isosurface extraction.
 This algorithm provides a good balance between performance and vertex count.
@@ -41,21 +22,13 @@ In contrast to the other algorithms, vertices may be repeated, so mesh size
 may be large and it will be difficult to extract topological/connectivity information.
 
 - `iso` (default: 0.0) specifies the iso level to use for surface extraction.
-- `eps` (default: 1e-3) is the tolerence around a voxel corner to ensure manifold mesh generation.
-- `reduceverts` (default: true) if true will merge vertices within a voxel to reduce mesh size by around 30% and with slight performance improvement.
 """
-struct MarchingCubes{T} <: AbstractMeshingAlgorithm
-    iso::T
-    eps::T
-    reduceverts::Bool
-    insidepositive::Bool
+Base.@kwdef struct MarchingCubes{T} <: AbstractMeshingAlgorithm
+    iso::T = 0.0
 end
 
 """
-    MarchingTetrahedra(iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
-    MarchingTetrahedra(;iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
-    MarchingTetrahedra(iso)
-    MarchingTetrahedra(iso,eps)
+    MarchingTetrahedra(;iso=0.0, eps=1e-3)
 
 Specifies the use of the Marching Tetrahedra algorithm for isosurface extraction.
 This algorithm has a roughly 2x performance penalty compared to Marching Cubes,
@@ -64,40 +37,15 @@ making this algorithm useful for topological analysis.
 
 - `iso` specifies the iso level to use for surface extraction.
 - `eps` is the tolerence around a voxel corner to ensure manifold mesh generation.
-- `reduceverts` (default: true) if false, vertices will not be unique and have repeated face references.
 """
-struct MarchingTetrahedra{T} <: AbstractMeshingAlgorithm
-    iso::T
-    eps::T
-    reduceverts::Bool
-    insidepositive::Bool
+Base.@kwdef struct MarchingTetrahedra{T,E} <: AbstractMeshingAlgorithm
+    iso::T = 0.0
+    eps::E = 1e-3
 end
 
 """
-    NaiveSurfaceNets(iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
-    NaiveSurfaceNets(;iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
-    NaiveSurfaceNets(iso)
-    NaiveSurfaceNets(iso,eps)
-
-Specifies the use of the Naive Surface Nets algorithm for isosurface extraction.
-This algorithm has a slight performance advantage in some cases over Marching Cubes.
-Each vertex is guaranteed to not be repeated (useful for topological analysis),
-however the algorithm does not guarantee accuracy and generates quad faces.
-
-- `iso` specifies the iso level to use for surface extraction.
-- `eps` is the tolerence around a voxel corner to ensure manifold mesh generation.
-- `reduceverts` reserved for future use.
-"""
-struct NaiveSurfaceNets{T} <: AbstractMeshingAlgorithm
-    iso::T
-    eps::T
-    reduceverts::Bool
-    insidepositive::Bool
-end
-
-"""
-    MarchingCubes(iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
-    MarchingCubes(;iso=0.0, eps=1e-3, reduceverts=true, insidepositive=false)
+    MarchingCubes(iso=0.0, eps=1e-3)
+    MarchingCubes(;iso=0.0, eps=1e-3)
     MarchingCubes(iso)
     MarchingCubes(iso,eps)
 
@@ -110,12 +58,11 @@ may be large and it will be difficult to extract topological/connectivity inform
 - `eps` (default: 1e-3) is the tolerence around a voxel corner to ensure manifold mesh generation.
 - `reduceverts` (default: true) if true will merge vertices within a voxel to reduce mesh size by around 30% and with slight performance improvement.
 """
-struct AdaptiveMarchingCubes{T} <: AbstractAdaptiveMeshingAlgorithm
+struct AdaptiveMarchingCubes{T,E} <: AbstractAdaptiveMeshingAlgorithm
     iso::T
     eps::T
     rtol::T
     atol::T
-    reduceverts::Bool
 end
 
 
@@ -124,13 +71,4 @@ struct AdaptiveMarchingTetrahedra{T} <: AbstractAdaptiveMeshingAlgorithm
     eps::T
     rtol::T
     atol::T
-    reduceverts::Bool
 end
-
-
-#
-# Helper functions
-#
-
-default_face_length(::Union{MarchingCubes,MarchingTetrahedra}) = 3
-default_face_length(::NaiveSurfaceNets) = 4
diff --git a/src/lut/mc.jl b/src/lut/mc.jl
index 570bfdd..86119c7 100644
--- a/src/lut/mc.jl
+++ b/src/lut/mc.jl
@@ -578,30 +578,30 @@ const _mc_eq_mapping = (0x01, 0x01, 0x02, 0x01, 0x03, 0x02, 0x04, 0x01, 0x02, 0x
                         0x15, 0x10, 0x10, 0x03, 0x0c, 0x14, 0x11, 0x01, 0x18, 0x0c, 0x0c,
                         0x11, 0x0c, 0x14, 0x11, 0x01, 0x0c, 0x11, 0x14, 0x01, 0x11, 0x01, 0x01)
 
-const _mc_connectivity = ((0x01, 0x02, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x04, 0x02, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x01, 0x04, 0x02, 0x04, 0x05, 0x02, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x03, 0x02, 0x04, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x01, 0x04, 0x02, 0x05, 0x06, 0x07, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x01, 0x03, 0x04, 0x01, 0x04, 0x05, 0x04, 0x03, 0x06, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x04, 0x01, 0x03, 0x04, 0x03, 0x05, 0x04, 0x05, 0x06, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x01, 0x04, 0x02, 0x01, 0x05, 0x04, 0x06, 0x02, 0x04, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x04, 0x06, 0x07, 0x06, 0x05, 0x08, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x01, 0x03, 0x04, 0x04, 0x03, 0x05, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x01, 0x03, 0x04, 0x01, 0x04, 0x05, 0x06, 0x04, 0x03, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x05, 0x07, 0x06, 0x05, 0x08, 0x07, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x03, 0x02, 0x04, 0x05, 0x06, 0x07, 0x05, 0x07, 0x08, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x01, 0x03, 0x04, 0x02, 0x05, 0x03, 0x06, 0x03, 0x07, 0x05, 0x07, 0x03),
-                            (0x01, 0x02, 0x03, 0x04, 0x01, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x0a, 0x0b, 0x0c, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x01, 0x04, 0x02, 0x01, 0x05, 0x04, 0x06, 0x04, 0x05, 0x07, 0x08, 0x09),
-                            (0x01, 0x02, 0x03, 0x01, 0x03, 0x04, 0x05, 0x06, 0x03, 0x06, 0x04, 0x03, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x01, 0x04, 0x02, 0x01, 0x05, 0x04, 0x04, 0x06, 0x02, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x04, 0x02, 0x01, 0x04, 0x05, 0x02, 0x04, 0x06, 0x05, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x01, 0x04, 0x02, 0x05, 0x06, 0x07, 0x06, 0x08, 0x07, 0x00, 0x00, 0x00),
-                            (0x01, 0x02, 0x03, 0x02, 0x04, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00))
+const _mc_connectivity = ((0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                        (0x04, 0x02, 0x01, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                        (0x04, 0x05, 0x06, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                        (0x01, 0x04, 0x02, 0x04, 0x05, 0x02, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                        (0x03, 0x02, 0x04, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                        (0x01, 0x04, 0x02, 0x05, 0x06, 0x07, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                        (0x01, 0x03, 0x04, 0x01, 0x04, 0x05, 0x04, 0x03, 0x06, 0x00, 0x00, 0x00),
+                        (0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                        (0x04, 0x01, 0x03, 0x04, 0x03, 0x05, 0x04, 0x05, 0x06, 0x00, 0x00, 0x00),
+                        (0x01, 0x04, 0x02, 0x01, 0x05, 0x04, 0x06, 0x02, 0x04, 0x00, 0x00, 0x00),
+                        (0x04, 0x05, 0x06, 0x04, 0x06, 0x07, 0x06, 0x05, 0x08, 0x00, 0x00, 0x00),
+                        (0x01, 0x03, 0x04, 0x04, 0x03, 0x05, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                        (0x01, 0x03, 0x04, 0x01, 0x04, 0x05, 0x06, 0x04, 0x03, 0x00, 0x00, 0x00),
+                        (0x04, 0x05, 0x06, 0x05, 0x07, 0x06, 0x05, 0x08, 0x07, 0x00, 0x00, 0x00),
+                        (0x03, 0x02, 0x04, 0x05, 0x06, 0x07, 0x05, 0x07, 0x08, 0x00, 0x00, 0x00),
+                        (0x01, 0x03, 0x04, 0x02, 0x05, 0x03, 0x06, 0x03, 0x07, 0x05, 0x07, 0x03),
+                        (0x04, 0x01, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00),
+                        (0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x0a, 0x0b, 0x0c, 0x00, 0x00, 0x00),
+                        (0x01, 0x04, 0x02, 0x01, 0x05, 0x04, 0x06, 0x04, 0x05, 0x07, 0x08, 0x09),
+                        (0x01, 0x03, 0x04, 0x05, 0x06, 0x03, 0x06, 0x04, 0x03, 0x00, 0x00, 0x00),
+                        (0x01, 0x04, 0x02, 0x01, 0x05, 0x04, 0x04, 0x06, 0x02, 0x00, 0x00, 0x00),
+                        (0x04, 0x02, 0x01, 0x04, 0x05, 0x02, 0x04, 0x06, 0x05, 0x00, 0x00, 0x00),
+                        (0x01, 0x04, 0x02, 0x05, 0x06, 0x07, 0x06, 0x08, 0x07, 0x00, 0x00, 0x00),
+                        (0x02, 0x04, 0x03, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00))
 
 const _mc_edge_list = ((0x01, 0x02), (0x02, 0x03), (0x03, 0x04), (0x04, 0x01),
                         (0x05, 0x06), (0x06, 0x07), (0x07, 0x08), (0x08, 0x05),
diff --git a/src/lut/mt.jl b/src/lut/mt.jl
index ae0e280..ad54018 100644
--- a/src/lut/mt.jl
+++ b/src/lut/mt.jl
@@ -20,20 +20,14 @@ Voxel corner and edge indexing conventions
   /
  X
 """
-# this gets vectorized so we want to ensure it is the
-# same type as out vertex
-@inline function voxCrnrPos(::Type{PT}) where {PT}
-    (PT(0, 0, 0),
-    PT(0, 1, 0),
-    PT(1, 1, 0),
-    PT(1, 0, 0),
-    PT(0, 0, 1),
-    PT(0, 1, 1),
-    PT(1, 1, 1),
-    PT(1, 0, 1))
-end
-
-const voxCrnrPosInt = voxCrnrPos(SVector{3,UInt8})
+const voxCrnrPos = ((0, 0, 0),
+                    (0, 1, 0),
+                    (1, 1, 0),
+                    (1, 0, 0),
+                    (0, 0, 1),
+                    (0, 1, 1),
+                    (1, 1, 1),
+                    (1, 0, 1))
 
 # the voxel IDs at either end of the tetrahedra edges, by edge ID
 const voxEdgeCrnrs = ((0x01, 0x02),
diff --git a/src/lut/sn.jl b/src/lut/sn.jl
deleted file mode 100644
index 610c62f..0000000
--- a/src/lut/sn.jl
+++ /dev/null
@@ -1,35 +0,0 @@
-const cube_edges = (0x00, 0x01, 0x00, 0x02, 0x00, 0x04, 0x01, 0x03, 0x01, 0x05, 0x02, 0x03,
-                    0x02, 0x06, 0x03, 0x07, 0x04, 0x05, 0x04, 0x06, 0x05, 0x07, 0x06, 0x07)
-
-const sn_edge_table = (0x0007, 0x0019, 0x001e, 0x0062, 0x0065, 0x007b, 0x007c, 0x00a8,
-                       0x00af, 0x00b1, 0x00b6, 0x00ca, 0x00cd, 0x00d3, 0x00d4, 0x0304,
-                       0x0303, 0x031d, 0x031a, 0x0366, 0x0361, 0x037f, 0x0378, 0x03ac,
-                       0x03ab, 0x03b5, 0x03b2, 0x03ce, 0x03c9, 0x03d7, 0x03d0, 0x0510,
-                       0x0517, 0x0509, 0x050e, 0x0572, 0x0575, 0x056b, 0x056c, 0x05b8,
-                       0x05bf, 0x05a1, 0x05a6, 0x05da, 0x05dd, 0x05c3, 0x05c4, 0x0614,
-                       0x0613, 0x060d, 0x060a, 0x0676, 0x0671, 0x066f, 0x0668, 0x06bc,
-                       0x06bb, 0x06a5, 0x06a2, 0x06de, 0x06d9, 0x06c7, 0x06c0, 0x0a40,
-                       0x0a47, 0x0a59, 0x0a5e, 0x0a22, 0x0a25, 0x0a3b, 0x0a3c, 0x0ae8,
-                       0x0aef, 0x0af1, 0x0af6, 0x0a8a, 0x0a8d, 0x0a93, 0x0a94, 0x0944,
-                       0x0943, 0x095d, 0x095a, 0x0926, 0x0921, 0x093f, 0x0938, 0x09ec,
-                       0x09eb, 0x09f5, 0x09f2, 0x098e, 0x0989, 0x0997, 0x0990, 0x0f50,
-                       0x0f57, 0x0f49, 0x0f4e, 0x0f32, 0x0f35, 0x0f2b, 0x0f2c, 0x0ff8,
-                       0x0fff, 0x0fe1, 0x0fe6, 0x0f9a, 0x0f9d, 0x0f83, 0x0f84, 0x0c54,
-                       0x0c53, 0x0c4d, 0x0c4a, 0x0c36, 0x0c31, 0x0c2f, 0x0c28, 0x0cfc,
-                       0x0cfb, 0x0ce5, 0x0ce2, 0x0c9e, 0x0c99, 0x0c87, 0x0c80, 0x0c80,
-                       0x0c87, 0x0c99, 0x0c9e, 0x0ce2, 0x0ce5, 0x0cfb, 0x0cfc, 0x0c28,
-                       0x0c2f, 0x0c31, 0x0c36, 0x0c4a, 0x0c4d, 0x0c53, 0x0c54, 0x0f84,
-                       0x0f83, 0x0f9d, 0x0f9a, 0x0fe6, 0x0fe1, 0x0fff, 0x0ff8, 0x0f2c,
-                       0x0f2b, 0x0f35, 0x0f32, 0x0f4e, 0x0f49, 0x0f57, 0x0f50, 0x0990,
-                       0x0997, 0x0989, 0x098e, 0x09f2, 0x09f5, 0x09eb, 0x09ec, 0x0938,
-                       0x093f, 0x0921, 0x0926, 0x095a, 0x095d, 0x0943, 0x0944, 0x0a94,
-                       0x0a93, 0x0a8d, 0x0a8a, 0x0af6, 0x0af1, 0x0aef, 0x0ae8, 0x0a3c,
-                       0x0a3b, 0x0a25, 0x0a22, 0x0a5e, 0x0a59, 0x0a47, 0x0a40, 0x06c0,
-                       0x06c7, 0x06d9, 0x06de, 0x06a2, 0x06a5, 0x06bb, 0x06bc, 0x0668,
-                       0x066f, 0x0671, 0x0676, 0x060a, 0x060d, 0x0613, 0x0614, 0x05c4,
-                       0x05c3, 0x05dd, 0x05da, 0x05a6, 0x05a1, 0x05bf, 0x05b8, 0x056c,
-                       0x056b, 0x0575, 0x0572, 0x050e, 0x0509, 0x0517, 0x0510, 0x03d0,
-                       0x03d7, 0x03c9, 0x03ce, 0x03b2, 0x03b5, 0x03ab, 0x03ac, 0x0378,
-                       0x037f, 0x0361, 0x0366, 0x031a, 0x031d, 0x0303, 0x0304, 0x00d4,
-                       0x00d3, 0x00cd, 0x00ca, 0x00b6, 0x00b1, 0x00af, 0x00a8, 0x007c,
-                       0x007b, 0x0065, 0x0062, 0x001e, 0x0019, 0x0007)
diff --git a/src/marching_cubes.jl b/src/marching_cubes.jl
index be4645e..4794048 100644
--- a/src/marching_cubes.jl
+++ b/src/marching_cubes.jl
@@ -2,28 +2,29 @@
 #Look up Table
 include("lut/mc.jl")
 
-function isosurface(sdf::AbstractArray{T, 3}, method::MarchingCubes, ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{3, Int};
-                    origin=VertType(-1,-1,-1), widths=VertType(2,2,2)) where {T, VertType, FaceType}
+function isosurface(sdf::AbstractArray{T,3}, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1) where {T}
     nx, ny, nz = size(sdf)
 
-    # we subtract one from the length along each axis because
-    # an NxNxN SDF has N-1 cells on each axis
-    s = VertType(widths[1]/(nx-1), widths[2]/(ny-1), widths[3]/(nz-1))
+    # find widest type
+    FT = promote_type(eltype(first(X)), eltype(first(Y)), eltype(first(Z)), eltype(T), typeof(method.iso))
 
-    # arrays for vertices and faces
-    vts = VertType[]
-    fcs = FaceType[]
+    vts = NTuple{3,float(FT)}[]
+    fcs = NTuple{3,Int}[]
+
+    xp = LinRange(first(X), last(X), nx)
+    yp = LinRange(first(Y), last(Y), ny)
+    zp = LinRange(first(Z), last(Z), nz)
 
     @inbounds for xi = 1:nx-1, yi = 1:ny-1, zi = 1:nz-1
 
-        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])
+        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
@@ -32,43 +33,41 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingCubes, ::Type{Vert
         # Cube is entirely in/out of the surface
         _no_triangles(cubeindex) && continue
 
-        points = mc_vert_points(xi,yi,zi,s,origin,VertType)
+        points = mc_vert_points(xi, yi, zi, xp, yp, zp)
 
-        # Create the triangle
-        _mc_unique_triangles!(method, points, iso_vals, vts, fcs, cubeindex, FaceType)
+        # process the voxel
+        process_mc_voxel!(vts, fcs, cubeindex, points, method.iso, iso_vals)
     end
-    vts,fcs
+    vts, fcs
 end
 
 
-function isosurface(f::Function, method::MarchingCubes, ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{3, Int};
-                    origin=VertType(-1,-1,-1), widths=VertType(2,2,2), samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {T <: Integer, VertType, FaceType}
+function isosurface(f::Function, method::MarchingCubes, X=-1:1, Y=-1:1, Z=-1:1; samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {T<:Integer}
 
     nx, ny, nz = samples[1], samples[2], samples[3]
 
-    # we subtract one from the length along each axis because
-    # an NxNxN SDF has N-1 cells on each axis
-    s = VertType(widths[1]/(nx-1), widths[2]/(ny-1), widths[3]/(nz-1))
+    # find widest type
+    FT = promote_type(eltype(first(X)), eltype(first(Y)), eltype(first(Z)), eltype(T), typeof(method.iso))
+
+    vts = NTuple{3,float(FT)}[]
+    fcs = NTuple{3,Int}[]
 
-    # arrays for vertices and faces
-    vts = VertType[]
-    fcs = FaceType[]
+    xp = LinRange(first(X), last(X), nx)
+    yp = LinRange(first(Y), last(Y), ny)
+    zp = LinRange(first(Z), last(Z), nz)
 
-    zv = zero(eltype(VertType))
-    iso_vals = (zv,zv,zv,zv,zv,zv,zv,zv)
     @inbounds for xi = 1:nx-1, yi = 1:ny-1, zi = 1:nz-1
 
-        points = mc_vert_points(xi,yi,zi,s,origin,VertType)
+        points = mc_vert_points(xi, yi, zi, xp, yp, zp)
 
         iso_vals = (f(points[1]),
-                    f(points[2]),
-                    f(points[3]),
-                    f(points[4]),
-                    f(points[5]),
-                    f(points[6]),
-                    f(points[7]),
-                    f(points[8]))
-        # TODO: Cache points
+            f(points[2]),
+            f(points[3]),
+            f(points[4]),
+            f(points[5]),
+            f(points[6]),
+            f(points[7]),
+            f(points[8]))
 
         #Determine the index into the edge table which
         #tells us which vertices are inside of the surface
@@ -77,83 +76,63 @@ function isosurface(f::Function, method::MarchingCubes, ::Type{VertType}=SVector
         # Cube is entirely in/out of the surface
         _no_triangles(cubeindex) && continue
 
-        # Create the triangle
-        _mc_unique_triangles!(method, points, iso_vals, vts, fcs, cubeindex, FaceType)
+        process_mc_voxel!(vts, fcs, cubeindex, points, method.iso, iso_vals)
     end
-    vts,fcs
+    vts, fcs
 end
 
-"""
-    _mc_unique_triangles!(vts, fcs, vertlist, cubeindex, FaceType)
-
-Create triangles by only adding unique vertices within the voxel.
-Each face may share a reference to a vertex with another face.
-"""
-function _mc_unique_triangles!(method, points, iso_vals, vts, fcs, cubeindex, ::Type{FaceType}) where {FaceType}
-    @inbounds begin
-        fct = length(vts)
 
-        find_vertices_interp!(vts, points, iso_vals, cubeindex, method.iso, method.eps)
-
-        offsets = _mc_connectivity[_mc_eq_mapping[cubeindex]]
-
-        # There is atleast one face so we can push it immediately
-        push!(fcs, FaceType(fct+offsets[3], fct+offsets[2], fct+offsets[1]))
-
-        for i in (4,7,10,13)
-            iszero(offsets[i]) && return
-            push!(fcs, FaceType(fct+offsets[i+2], fct+offsets[i+1], fct+offsets[i]))
-        end
-    end
-end
+function process_mc_voxel!(vts, fcs, cubeindex, points, iso, iso_vals)
 
+    fct = length(vts)
 
-"""
-    find_vertices_interp(points, iso_vals, cubeindex, iso, eps)
-
-Find the vertices where the surface intersects the cube
-"""
-function find_vertices_interp!(vts, points, iso_vals, cubeindex, iso, eps)
     @inbounds begin
+        # Add the vertices
         vert_to_add = _mc_verts[cubeindex]
         for i = 1:12
             vt = vert_to_add[i]
             iszero(vt) && break
             ed = _mc_edge_list[vt]
-            push!(vts, vertex_interp(iso,points[ed[1]],points[ed[2]],iso_vals[ed[1]],iso_vals[ed[2]], eps))
+            push!(vts, vertex_interp(iso, points[ed[1]], points[ed[2]], iso_vals[ed[1]], iso_vals[ed[2]]))
+        end
+
+        # Add the faces
+        offsets = _mc_connectivity[_mc_eq_mapping[cubeindex]]
+
+        # There is atleast one face so we can push it immediately
+        push!(fcs, (fct + 3, fct + 2, fct + 1))
+
+        for i in (1, 4, 7, 10)
+            iszero(offsets[i]) && return
+            push!(fcs, (fct + offsets[i+2], fct + offsets[i+1], fct + offsets[i]))
         end
     end
 end
 
-"""
-    vertex_interp(iso, p1, p2, valp1, valp2, eps = 0.00001)
+
+"""    vertex_interp(iso, p1, p2, valp1, valp2)
 
 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)
-
-    abs(iso - valp1) < eps && return p1
-    abs(iso - valp2) < eps && return p2
-    abs(valp1-valp2) < eps && return p1
+function vertex_interp(iso, p1, p2, valp1, valp2)
     mu = (iso - valp1) / (valp2 - valp1)
-    p = p1 + mu * (p2 - p1)
-
+    p = p1 .+ mu .* (p2 .- p1)
     return p
 end
 
 """
-    mc_vert_points(xi,yi,zi, scale, origin, ::Type{VertType})
+    mc_vert_points(xi,yi,zi,xp,yp,zp)
 
 Returns a tuple of 8 points corresponding to each corner of a cube
 """
-@inline function mc_vert_points(xi,yi,zi, scale, origin, ::Type{VertType}) where VertType
-    (VertType(xi-1,yi-1,zi-1) .* scale .+ origin,
-     VertType(xi  ,yi-1,zi-1) .* scale .+ origin,
-     VertType(xi  ,yi  ,zi-1) .* scale .+ origin,
-     VertType(xi-1,yi  ,zi-1) .* scale .+ origin,
-     VertType(xi-1,yi-1,zi  ) .* scale .+ origin,
-     VertType(xi  ,yi-1,zi  ) .* scale .+ origin,
-     VertType(xi  ,yi  ,zi  ) .* scale .+ origin,
-     VertType(xi-1,yi  ,zi  ) .* scale .+ origin)
+function mc_vert_points(xi, yi, zi, xp, yp, zp)
+    ((xp[xi], yp[yi], zp[zi]),
+        (xp[xi+1], yp[yi], zp[zi]),
+        (xp[xi+1], yp[yi+1], zp[zi]),
+        (xp[xi], yp[yi+1], zp[zi]),
+        (xp[xi], yp[yi], zp[zi+1]),
+        (xp[xi+1], yp[yi], zp[zi+1]),
+        (xp[xi+1], yp[yi+1], zp[zi+1]),
+        (xp[xi], yp[yi+1], zp[zi+1]))
 end
diff --git a/src/marching_tetrahedra.jl b/src/marching_tetrahedra.jl
index aacccc5..faffc4a 100644
--- a/src/marching_tetrahedra.jl
+++ b/src/marching_tetrahedra.jl
@@ -8,7 +8,7 @@ include("lut/mt.jl")
 
 Determines which case in the triangle table we are dealing with
 """
-@inline function tetIx(tIx, cubeindex)
+function tetIx(tIx, cubeindex)
     @inbounds v1 = subTetsMask[tIx][1]
     @inbounds v2 = subTetsMask[tIx][2]
     ix = 0x01 & cubeindex | (0x40 & cubeindex) >> 0x03
@@ -26,8 +26,8 @@ two edges get the same index) and unique (every edge gets the same ID
 regardless of which of its neighboring voxels is asking for it) in order
 for vertex sharing to be implemented properly.
 """
-@inline function vertId(e, x, y, z, nx, ny)
-    dx, dy, dz = voxCrnrPosInt[voxEdgeCrnrs[e][1]]
+function vertId(e, x, y, z, nx, ny)
+    dx, dy, dz = voxCrnrPos[voxEdgeCrnrs[e][1]]
     voxEdgeDir[e]+7*(x-0x01+dx+nx*(y-0x01+dy+ny*(z-0x01+dz)))
 end
 
@@ -39,7 +39,7 @@ occurs.
 eps represents the "bump" factor to keep vertices away from voxel
 corners (thereby preventing degeneracies).
 """
-@inline function vertPos(e, x, y, z, scale, origin, vals::V, iso, eps, ::Type{VertType}) where {V, VertType}
+function vertPos(e, x, y, z, xp, yp, zp, vals::V, iso, eps) where {V}
     T = eltype(vals)
 
     ixs     = voxEdgeCrnrs[e]
@@ -47,10 +47,11 @@ corners (thereby preventing degeneracies).
     tgtVal  = vals[ixs[2]]
     a       = min(max((iso-srcVal)/(tgtVal-srcVal), eps), one(T)-eps)
     b       = one(T)-a
-    c1 = voxCrnrPos(VertType)[ixs[1]]
-    c2 = voxCrnrPos(VertType)[ixs[2]]
+    c1 = voxCrnrPos[ixs[1]]
+    c2 = voxCrnrPos[ixs[2]]
 
-    ((VertType(x,y,z) + c1 .* b + c2.* a) .- 1) .* scale .+ origin
+    d = (xp[x+1], yp[y+1], zp[z+1]) .- (xp[x], yp[y], zp[z])
+    (xp[x],yp[y],zp[z]) .+ (c1 .* b .+ c2 .* a) .* d
 end
 
 """
@@ -63,28 +64,21 @@ end
 Gets the vertex ID, adding it to the vertex dictionary if not already
 present.
 """
-@inline function getVertId(e, x, y, z, nx, ny,
+function getVertId(e, x, y, z, nx, ny,
                            vals, iso::Real,
-                           scale, origin,
+                           xp, yp, zp,
                            vts::Dict,
                            vtsAry::Vector,
-                           eps::Real,
-                           reduceverts::Bool)
+                           eps::Real)
 
-    VertType = eltype(vtsAry)
-    if reduceverts
-        vId = vertId(e, x, y, z, nx, ny)
-        haskey(vts, vId) && return vts[vId]
-    end
+    vId = vertId(e, x, y, z, nx, ny)
+    haskey(vts, vId) && return vts[vId]
 
     # calculate vert position
-    v = vertPos(e, x, y, z, scale, origin, vals, iso, eps, VertType)
+    v = vertPos(e, x, y, z, xp, yp, zp, vals, iso, eps)
     push!(vtsAry, v)
 
-    # if deduplicting, push to dict
-    if reduceverts
-        vts[vId] = length(vtsAry)
-    end
+    vts[vId] = length(vtsAry)
 
     return length(vtsAry)
 end
@@ -94,7 +88,7 @@ end
 Given a sub-tetrahedron case and a tetrahedron edge ID, determines the
 corresponding voxel edge ID.
 """
-@inline function voxEdgeId(subTetIx, tetEdgeIx)
+function voxEdgeId(subTetIx, tetEdgeIx)
     srcVoxCrnr = subTets[subTetIx][tetEdgeCrnrs[tetEdgeIx][1]]
     tgtVoxCrnr = subTets[subTetIx][tetEdgeCrnrs[tetEdgeIx][2]]
     return voxEdgeIx[srcVoxCrnr][tgtVoxCrnr]
@@ -108,11 +102,10 @@ end
 Processes a voxel, adding any new vertices and faces to the given
 containers as necessary.
 """
-function procVox(vals, iso::Real, x, y, z, nx, ny, scale, origin,
+function procVox(vals, iso::Real, x, y, z, nx, ny, xp, yp, zp,
                  vts::Dict, vtsAry::Vector, fcs::Vector,
-                 eps::Real, cubeindex, reduceverts)
-    VertType = eltype(vtsAry)
-    FaceType = eltype(fcs)
+                 eps::Real, cubeindex)
+
     # check each sub-tetrahedron in the voxel
     @inbounds for i = 1:6
         tIx = tetIx(i, cubeindex)
@@ -121,31 +114,33 @@ function procVox(vals, iso::Real, x, y, z, nx, ny, scale, origin,
         e = tetTri[tIx]
 
         # add the face to the list
-        push!(fcs, FaceType(
-                    getVertId(voxEdgeId(i, e[1]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts),
-                    getVertId(voxEdgeId(i, e[2]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts),
-                    getVertId(voxEdgeId(i, e[3]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts)))
+        push!(fcs, (getVertId(voxEdgeId(i, e[1]), x, y, z, nx, ny, vals, iso, xp, yp, zp, vts, vtsAry, eps),
+                    getVertId(voxEdgeId(i, e[2]), x, y, z, nx, ny, vals, iso, xp, yp, zp, vts, vtsAry, eps),
+                    getVertId(voxEdgeId(i, e[3]), x, y, z, nx, ny, vals, iso, xp, yp, zp, vts, vtsAry, eps)))
 
         # bail if there are no more faces
         iszero(e[4]) && continue
-        push!(fcs, FaceType(
-                    getVertId(voxEdgeId(i, e[4]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts),
-                    getVertId(voxEdgeId(i, e[5]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts),
-                    getVertId(voxEdgeId(i, e[6]), x, y, z, nx, ny, vals, iso, scale, origin, vts, vtsAry, eps, reduceverts)))
+        push!(fcs, (getVertId(voxEdgeId(i, e[4]), x, y, z, nx, ny, vals, iso, xp, yp, zp, vts, vtsAry, eps),
+                    getVertId(voxEdgeId(i, e[5]), x, y, z, nx, ny, vals, iso, xp, yp, zp, vts, vtsAry, eps),
+                    getVertId(voxEdgeId(i, e[6]), x, y, z, nx, ny, vals, iso, xp, yp, zp, vts, vtsAry, eps)))
     end
 end
 
-function isosurface(sdf::AbstractArray{T, 3}, method::MarchingTetrahedra, ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{3, Int};
-                    origin=VertType(-1,-1,-1), widths=VertType(2,2,2)) where {T, VertType, FaceType}
+function isosurface(sdf::AbstractArray{T, 3}, method::MarchingTetrahedra, X=-1:1, Y=-1:1, Z=-1:1) where {T}
+
+    FT = promote_type(eltype(first(X)), eltype(first(Y)), eltype(first(Z)), eltype(T), typeof(method.iso), typeof(method.eps))
 
     vts    = Dict{Int, Int}()
-    fcs    = FaceType[]
-    vtsAry = VertType[]
+    fcs = NTuple{3,Int}[]
+    vtsAry = NTuple{3,float(FT)}[]
 
     # process each voxel
-    scale = widths ./ VertType(size(sdf) .- 1)
     nx::Int, ny::Int, nz::Int = size(sdf)
 
+    xp = LinRange(first(X), last(X), nx)
+    yp = LinRange(first(Y), last(Y), ny)
+    zp = LinRange(first(Z), last(Z), nz)
+
     @inbounds for i = 1:nx-1, j = 1:ny-1, k = 1:nz-1
 
         vals = (sdf[i,  j,  k  ],
@@ -161,36 +156,30 @@ function isosurface(sdf::AbstractArray{T, 3}, method::MarchingTetrahedra, ::Type
 
         _no_triangles(cubeindex) && continue
 
-        procVox(vals, method.iso, i, j, k, nx, ny, scale, origin, vts, vtsAry, fcs, method.eps, cubeindex, method.reduceverts)
+        procVox(vals, method.iso, i, j, k, nx, ny, xp, yp, zp, vts, vtsAry, fcs, method.eps, cubeindex)
     end
 
     vtsAry,fcs
 end
 
-function isosurface(f::Function, method::MarchingTetrahedra,
-                    ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{3, Int};
-                    origin=VertType(-1,-1,-1), widths=VertType(2,2,2),
-                    samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {T, VertType, FaceType}
+function isosurface(f::F, method::MarchingTetrahedra, X=-1:1, Y=-1:1, Z=-1:1;
+                    samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {F, T}
+
+    FT = promote_type(eltype(first(X)), eltype(first(Y)), eltype(first(Z)), eltype(T), typeof(method.iso), typeof(method.eps))
 
     vts    = Dict{Int, Int}()
-    fcs    = FaceType[]
-    vtsAry = VertType[]
+    fcs    = NTuple{3,Int}[]
+    vtsAry = NTuple{3,float(FT)}[]
 
     # process each voxel
-    scale = widths ./ (VertType(samples...) .-1)
-    nx::Int, ny::Int, nz::Int = samples
-    zv = zero(eltype(VertType))
-    vals = (zv,zv,zv,zv,zv,zv,zv,zv)
+    nx, ny, nz = samples[1], samples[2], samples[3]
+    xp = LinRange(first(X), last(X), nx)
+    yp = LinRange(first(Y), last(Y), ny)
+    zp = LinRange(first(Z), last(Z), nz)
 
     @inbounds for i = 1:nx-1, j = 1:ny-1, k = 1:nz-1
-        points = (VertType(i-1,j-1,k-1).* scale .+ origin,
-                  VertType(i-1,j  ,k-1).* scale .+ origin,
-                  VertType(i  ,j  ,k-1).* scale .+ origin,
-                  VertType(i  ,j-1,k-1).* scale .+ origin,
-                  VertType(i-1,j-1,k  ).* scale .+ origin,
-                  VertType(i-1,j  ,k  ).* scale .+ origin,
-                  VertType(i  ,j  ,k  ).* scale .+ origin,
-                  VertType(i  ,j-1,k  ).* scale .+ origin)
+        points = mt_vert_points(i, j, k, xp, yp, zp)
+
         vals = (f(points[1]),
                 f(points[2]),
                 f(points[3]),
@@ -204,8 +193,19 @@ function isosurface(f::Function, method::MarchingTetrahedra,
 
         _no_triangles(cubeindex) && continue
 
-        procVox(vals, method.iso, i, j, k, nx, ny, scale, origin, vts, vtsAry, fcs, method.eps, cubeindex, method.reduceverts)
+        procVox(vals, method.iso, i, j, k, nx, ny, xp, yp, zp, vts, vtsAry, fcs, method.eps, cubeindex)
     end
 
     vtsAry,fcs
 end
+
+function mt_vert_points(i, j, k, xp, yp, zp)
+    ((xp[i], yp[j], zp[k]),
+    (xp[i], yp[j+1], zp[k]),
+    (xp[i+1], yp[j+1], zp[k]),
+    (xp[i+1], yp[j], zp[k]),
+    (xp[i], yp[j], zp[k+1]),
+    (xp[i], yp[j+1], zp[k+1]),
+    (xp[i+1], yp[j+1], zp[k+1]),
+    (xp[i+1], yp[j], zp[k+1]))
+end
\ No newline at end of file
diff --git a/src/surface_nets.jl b/src/surface_nets.jl
deleted file mode 100644
index 5f618d0..0000000
--- a/src/surface_nets.jl
+++ /dev/null
@@ -1,270 +0,0 @@
-#
-# SurfaceNets in Julia
-#
-# Ported from the Javascript implementation by Mikola Lysenko (C) 2012
-# https://github.com/mikolalysenko/mikolalysenko.github.com/blob/master/Isosurface/js/surfacenets.js
-# MIT License
-#
-# Based on: S.F. Gibson, "Constrained Elastic Surface Nets". (1998) MERL Tech Report.
-#
-
-#Look up Table
-include("lut/sn.jl")
-
-function isosurface(sdf::AbstractArray{T, 3}, method::NaiveSurfaceNets, ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{4, Int};
-                    origin=VertType(-1,-1,-1), widths=VertType(2,2,2)) where {T, VertType, FaceType}
-
-    scale = widths ./ VertType(size(sdf) .- 1)  # subtract 1 because an SDF with N points per side has N-1 cells
-
-    dims = size(sdf)
-
-    vertices = VertType[]
-    faces = FaceType[]
-
-    n = 0
-    R = Array{Int}([1, (dims[1]+1), (dims[1]+1)*(dims[2]+1)])
-    buf_no = 1
-
-    buffer = fill(zero(Int),R[3]*2)
-
-    #March over the voxel grid
-    zi = 0
-    @inbounds while zi<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))
-
-        for yi = 0:dims[2]-2
-            for xi = 0:dims[1]-2
-
-                inds = (xi,yi,zi)
-
-                # Read in 8 field values around this vertex and store them in an array
-                # Also calculate 8-bit mask, like in marching cubes, so we can speed up sign checks later
-                @inbounds grid = (sdf[xi+1,yi+1,zi+1],
-                                  sdf[xi+2,yi+1,zi+1],
-                                  sdf[xi+1,yi+2,zi+1],
-                                  sdf[xi+2,yi+2,zi+1],
-                                  sdf[xi+1,yi+1,zi+2],
-                                  sdf[xi+2,yi+1,zi+2],
-                                  sdf[xi+1,yi+2,zi+2],
-                                  sdf[xi+2,yi+2,zi+2])
-
-                mask = _get_cubeindex(grid, method.iso)
-
-                # Check for early termination if cell does not intersect boundary
-                if _no_triangles(mask)
-                    m += 1
-                    continue
-                end
-
-                # iso level correction
-                grid = grid .- method.iso
-
-                #Sum up edge intersections
-                edge_mask = sn_edge_table[mask]
-
-                _sn_add_verts!(inds, vertices, grid, edge_mask, buffer, m, scale, origin, method.eps, VertType)
-
-                #Now we need to add faces together, to do this we just loop over 3 basis components
-                _sn_add_faces!(inds, faces, edge_mask, mask, buffer, m, R, FaceType)
-
-                m += 1
-            end
-            m += 2
-        end
-        zi += 1
-        buf_no = xor(buf_no,1)
-        R[3]=-R[3]
-    end
-
-    vertices, faces
-end
-
-function isosurface(f::Function, method::NaiveSurfaceNets,
-                    ::Type{VertType}=SVector{3,Float64}, ::Type{FaceType}=SVector{4, Int};
-                    origin=VertType(-1,-1,-1), widths=VertType(2,2,2),
-                    samples::NTuple{3,T}=_DEFAULT_SAMPLES) where {T <: Integer, VertType, FaceType}
-
-    scale = widths ./ VertType(samples .- 1)  # subtract 1 because an SDF with N points per side has N-1 cells
-
-    vertices = VertType[]
-    faces = FaceType[]
-
-
-    R = Array{Int}([1, (samples[1]+1), (samples[1]+1)*(samples[2]+1)])
-    buf_no = 1
-
-    buffer = fill(zero(Int),R[3]*2)
-
-    zv = zero(eltype(VertType))
-    zvt = zero(VertType)
-    grid = (zv,zv,zv,zv,zv,zv,zv,zv)
-    points = (zvt,zvt,zvt,zvt,zvt,zvt,zvt,zvt)
-
-    #March over the voxel grid
-    zi = 0
-    @inbounds while zi<samples[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 + (samples[1]+1) * (1 + buf_no * (samples[2]+1))
-
-        for yi = 0:samples[2]-2
-            for xi = 0:samples[1]-2
-
-                inds = (xi,yi,zi)
-
-                if xi == 0
-                    points = (VertType(xi,yi,zi).* scale .+ origin,
-                              VertType(xi+1,yi,zi).* scale .+ origin,
-                              VertType(xi,yi+1,zi).* scale .+ origin,
-                              VertType(xi+1,yi+1,zi).* scale .+ origin,
-                              VertType(xi,yi,zi+1).* scale .+ origin,
-                              VertType(xi+1,yi,zi+1).* scale .+ origin,
-                              VertType(xi,yi+1,zi+1).* scale .+ origin,
-                              VertType(xi+1,yi+1,zi+1).* scale .+ origin)
-                    grid = (f(points[1]),
-                            f(points[2]),
-                            f(points[3]),
-                            f(points[4]),
-                            f(points[5]),
-                            f(points[6]),
-                            f(points[7]),
-                            f(points[8]))
-                else
-                    points = (points[2],
-                              VertType(xi+1,yi,zi).* scale .+ origin,
-                              points[4],
-                              VertType(xi+1,yi+1,zi).* scale .+ origin,
-                              points[6],
-                              VertType(xi+1,yi,zi+1).* scale .+ origin,
-                              points[8],
-                              VertType(xi+1,yi+1,zi+1).* scale .+ origin)
-                    grid = (grid[2],
-                            f(points[2]),
-                            grid[4],
-                            f(points[4]),
-                            grid[6],
-                            f(points[6]),
-                            grid[8],
-                            f(points[8]))
-                end
-
-                # Also calculate 8-bit mask, like in marching cubes, so we can speed up sign checks later
-                mask = _get_cubeindex(grid, method.iso)
-
-                # Check for early termination if cell does not intersect boundary
-                if _no_triangles(mask)
-                    m += 1
-                    continue
-                end
-
-                # iso level correction
-                grid = grid .- method.iso
-
-                #Sum up edge intersections
-                edge_mask = sn_edge_table[mask]
-
-                # add vertices
-                _sn_add_verts!(inds, vertices, grid, edge_mask, buffer, m, scale, origin, method.eps, VertType)
-
-                #Now we need to add faces together, to do this we just loop over 3 basis components
-                _sn_add_faces!(inds, faces, edge_mask, mask, buffer, m, R, FaceType)
-
-                m += 1
-            end
-            m += 2
-        end
-        zi += 1
-        buf_no = xor(buf_no,1)
-        R[3]=-R[3]
-    end
-
-    vertices, faces
-end
-
-@inline function _sn_add_verts!(inds, vertices, grid, edge_mask, buffer, m, scale, origin, eps, ::Type{VertType}) where {VertType}
-    v = zero(VertType)
-    T = eltype(VertType)
-    e_count = 0
-
-    #For every edge of the cube...
-    @inbounds for i=0x00:0x0b
-
-        #Use edge mask to check if it is crossed
-        iszero(edge_mask & (0x0001<<i)) && continue
-
-        #If it did, increment number of edge crossings
-        e_count += 1
-
-        #Now find the point of intersection
-        e0 = cube_edges[(i<<0x01)+1]       #Unpack vertices
-        e1 = cube_edges[(i<<0x01)+2]
-        g0 = grid[e0+1]                 #Unpack grid values
-        g1 = grid[e1+1]
-        t  = g0 - g1                 #Compute point of intersection
-
-        abs(t) <= eps && continue
-
-        t = g0 / t
-
-        #Interpolate vertices and add up intersections (this can be done without multiplying)
-        a1, a2, a3 = 0x01 & e0, 0x02 & e0, 0x04 & e0
-        b1, b2, b3 = 0x01 & e1, 0x02 & e1, 0x04 & e1
-        xj = T(!iszero(a1))
-        yj = T(!iszero(a2))
-        zj = T(!iszero(a3))
-        a1 != b1 && (xj += !iszero(a1) ? -t : t)
-        a2 != b2 && (yj += !iszero(a2) ? -t : t)
-        a3 != b3 && (zj += !iszero(a3) ? -t : t)
-        v = v .+ VertType(xj,yj,zj)
-
-    end # edge check
-
-    #Now we just average the edge intersections and add them to coordinate
-    s = one(T) / e_count
-    v = (VertType(inds...)  .+ s .* v) .* scale .+ origin
-
-    #Add vertex to buffer, store pointer to vertex index in buffer
-    buffer[m+1] = length(vertices)
-    push!(vertices, v)
-end
-
-function _sn_add_faces!(inds, faces, edge_mask, mask, buffer, m, R, ::Type{FaceType}) where {FaceType}
-    for i = 0x00:0x02
-        #The first three entries of the edge_mask count the crossings along the edge
-        iszero(edge_mask & (0x0001<<i)) && continue
-
-        # i = axes we are point along.  iu, iv = orthogonal axes
-        iu = (i+0x01)%0x03
-        iv = (i+0x02)%0x03
-
-        #If we are on a boundary, skip it
-        iszero(inds[iu+1]) || iszero(inds[iv+1]) && continue
-
-        #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 !iszero(mask & 0x01)
-            if length(FaceType) == 4
-                push!(faces,FaceType(buffer[m+1]+1, buffer[m-du+1]+1, buffer[m-du-dv+1]+1, buffer[m-dv+1]+1))
-            elseif length(FaceType) == 3
-                push!(faces,FaceType(buffer[m+1]+1, buffer[m-du+1]+1, buffer[m-du-dv+1]+1))
-                push!(faces,FaceType(buffer[m-du-dv+1]+1, buffer[m-dv+1]+1, buffer[m+1]+1))
-            end
-        else
-            if length(FaceType) == 4
-                push!(faces,FaceType(buffer[m+1]+1, buffer[m-dv+1]+1, buffer[m-du-dv+1]+1, buffer[m-du+1]+1))
-            elseif length(FaceType) == 3
-                push!(faces,FaceType(buffer[m+1]+1, buffer[m-dv+1]+1, buffer[m-du-dv+1]+1))
-                push!(faces,FaceType(buffer[m-du-dv+1]+1, buffer[m-du+1]+1, buffer[m+1]+1))
-            end
-        end
-    end
-end
diff --git a/test/runtests.jl b/test/runtests.jl
index ebe0d91..6073c62 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,26 +1,14 @@
 using Meshing
 using Test
 using ForwardDiff
-using StaticArrays
 using Statistics: mean
 using LinearAlgebra: dot, norm
 using Random
 
-const algos = (MarchingCubes, MarchingTetrahedra, NaiveSurfaceNets)
+const algos = (MarchingCubes, MarchingTetrahedra)
 
 sphere_function(v) = sqrt(sum(dot(v, v))) - 1
 torus_function(v) = (sqrt(v[1]^2 + v[2]^2) - 0.5)^2 + v[3]^2 - 0.25
-@testset "meshing algorithms" begin
-    for algo in (MarchingCubes, MarchingTetrahedra, NaiveSurfaceNets)
-        @test algo(5) == algo{Float64}(5.0, 0.001, true, false)
-        @test algo(0x00, 0x00) == algo{UInt8}(0x00, 0x00, true, false)
-        @test algo{UInt16}(5, 0x00) == algo{UInt16}(0x0005, 0x0000, true, false)
-        @test algo{Float32}(1) ==  algo{Float32}(1.0f0, 0.001f0, true, false)
-        @test algo{Float32}(1.0, 0x00) == algo{Float32}(1.0f0, 0.0f0, true, false)
-        @test algo{Float32}(iso=1, eps=0x00, insidepositive=true, reduceverts=false) ==
-            algo{Float32}(1.0f0, 0.0f0, false, true)
-    end
-end
 
 @testset "_get_cubeindex" begin
     iso_vals1 = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8]
@@ -36,62 +24,38 @@ end
     @test Meshing._get_cubeindex(iso_vals3, iso3) == 0xe0
 end
 
+@testset "vertex interpolation" begin
+    @test Meshing.vertex_interp(0, (0,0,0), (0,1,0), -1, 1) == (0,0.5,0)
+    @test Meshing.vertex_interp(-1, (0,0,0), (0,1,0), -1, 1) == (0,0,0)
+    @test Meshing.vertex_interp(1, (0,0,0), (0,1,0), -1, 1) == (0,1,0)
+end
+
 @testset "respect origin" begin
     # verify that when we construct a mesh, that mesh:
     # a) respects the origin of the SDF
     # b) is correctly spaced so that symmetric functions create symmetric meshes
     f = x -> norm(x)
-    origin = SVector{3, Float64}(-1, -1, -1)
-    widths = SVector{3, Float64}(2, 2, 2)
     resolution = 0.1
 
-    for algorithm in (MarchingCubes(0.5),
-                      MarchingTetrahedra(0.5),
-                      MarchingCubes(iso=0.5, reduceverts=false),
-                      MarchingTetrahedra(iso=0.5, reduceverts=false))
+    for algorithm in (MarchingCubes(iso=0.5),
+                      MarchingTetrahedra(iso=0.5))
         @testset "$(typeof(algorithm))" begin
             # Extract isosurface using a function
-            points, faces = isosurface(f, algorithm, origin=origin, widths=widths, samples=(50, 50, 50))
+            points, faces = isosurface(f, algorithm, samples=(50, 50, 50))
 
             # should be centered on the origin
-            @test mean(points) ≈ [0, 0, 0] atol=0.15*resolution
+            @test mean(collect.(points)) ≈ [0, 0, 0] atol=0.15*resolution
             # and should be symmetric about the origin
-            #@test maximum(points) ≈ [0.5, 0.5, 0.5]
-            #@test minimum(points) ≈ [-0.5, -0.5, -0.5]
+
+            for i in 1:3
+                @test maximum(map(p -> p[i], collect.(points))) ≈ 0.5 atol=0.001
+                @test minimum(map(p -> p[i], collect.(points))) ≈ -0.5 atol=0.001
+            end
         end
     end
 
-    #TODO: SurfaceNets functional variant is bugged
-    # Naive Surface Nets has no accuracy guarantee, and is a weighted sum
-    # so a larger tolerance is needed for this one. In addition,
-    # quad -> triangle conversion is not functioning correctly
-    # see: https://github.com/JuliaGeometry/GeometryTypes.jl/issues/169
-    points, faces = isosurface(f, NaiveSurfaceNets(0.5), origin=origin, widths=widths, samples=(21, 21, 21))
-
-    # should be centered on the origin
-    #@test mean(points) ≈ [0, 0, 0] atol=0.15*resolution
-    # and should be symmetric about the origin
-    #@test maximum(points) ≈ [0.5, 0.5, 0.5] atol=0.2
-    #@test minimum(points) ≈ [-0.5, -0.5, -0.5] atol=0.2
 end
-@testset "surface nets" begin
-
 
-    """
-    Test the isosurface function with NaiveSurfaceNets and a function.
-
-    This code tests the `isosurface` function with the `NaiveSurfaceNets` algorithm and two different level set functions: `sphere_function` and `torus_function`. It checks that the number of vertices and faces in the resulting meshes match expected values.
-    """
-    # Test the isosurface function with NaiveSurfaceNets and a function
-    points_sphere, faces_sphere = isosurface(sphere_function, NaiveSurfaceNets(), origin=SVector(-1, -1, -1), widths=SVector(2, 2, 2), samples=(21, 21, 21))
-    points_torus, faces_torus = isosurface(torus_function, NaiveSurfaceNets(), origin=SVector(-2, -2, -2), widths=SVector(4, 4, 4), samples=(81, 81, 81))
-
-    # Add assertions to check the output
-    @test length(points_sphere) == 1832
-    @test length(faces_sphere) == 1830
-    @test length(points_torus) == 5532
-    @test length(faces_torus) == 5532
-end
 
 @testset "noisy spheres" begin
     # Produce a level set function that is a noisy version of the distance from
@@ -109,23 +73,17 @@ end
     #
     lambda = N-2*sigma # isovalue
 
-    points, faces = isosurface(distance, MarchingTetrahedra(lambda))
+    points, faces = isosurface(distance, MarchingTetrahedra(iso=lambda))
 
     @test length(points) == 3466
     @test length(faces) == 6928
 end
 
-@testset "vertex interpolation" begin
-    @test Meshing.vertex_interp(0, SVector(0,0,0), SVector(0,1,0), -1, 1) == SVector(0,0.5,0)
-    @test Meshing.vertex_interp(-1, SVector(0,0,0), SVector(0,1,0), -1, 1) == SVector(0,0,0)
-    @test Meshing.vertex_interp(1, SVector(0,0,0), SVector(0,1,0), -1, 1) == SVector(0,1,0)
-end
-
 @testset "marching cubes" begin
     algo = MarchingCubes()
 
     # Extract isosurfaces using MarchingCubes
-    points_mf, faces_mf = isosurface(sphere_function, algo, origin=SVector(-1, -1, -1), widths=SVector(2, 2, 2), samples=(21, 21, 21))
+    points_mf, faces_mf = isosurface(sphere_function, algo, samples=(21, 21, 21))
 
     # Test the number of vertices and faces
     @test length(points_mf) == 7320
@@ -133,27 +91,20 @@ end
 end
 @testset "marching tetrahedra" begin
     a1 = MarchingTetrahedra()
-    a2 = MarchingTetrahedra(reduceverts=false)
 
     # Extract isosurfaces using MarchingTetrahedra
-    points_mt1, faces_mt1 = isosurface(v -> sqrt(sum(dot(v, v))) - 1, a1, origin=SVector(-1, -1, -1), widths=SVector(2, 2, 2), samples=(21, 21, 21))
-    points_mt2, faces_mt2 = isosurface(v -> sqrt(sum(dot(v, v))) - 1, a2, origin=SVector(-1, -1, -1), widths=SVector(2, 2, 2), samples=(21, 21, 21))
+    points_mt1, faces_mt1 = isosurface(v -> sqrt(sum(dot(v, v))) - 1, a1, samples=(21, 21, 21))
 
     # Test the number of vertices and faces
     @test length(points_mt1) == 5582
-    @test length(points_mt2) == 33480
     @test length(faces_mt1) == 11160
-    @test length(faces_mt2) == length(faces_mt1)
 
     # When zero set is not completely within the box
-    points_mt1_partial, faces_mt1_partial = isosurface(v -> v[3] - 50, a1, origin=SVector(0, 0, 40), widths=SVector(50, 50, 60), samples=(2, 2, 2))
-    points_mt2_partial, faces_mt2_partial = isosurface(v -> v[3] - 50, a2, origin=SVector(0, 0, 40), widths=SVector(50, 50, 60), samples=(2, 2, 2))
+    points_mt1_partial, faces_mt1_partial = isosurface(v -> v[3] - 50, a1, 0:50, 0:50, 40:60, samples=(2, 2, 2))
 
     # Test the number of vertices and faces
     @test length(points_mt1_partial) == 9
-    @test length(points_mt2_partial) == 24
     @test length(faces_mt1_partial) == 8
-    @test length(faces_mt2_partial) == length(faces_mt1_partial)
 end
 
 @testset "function/array" begin
@@ -161,11 +112,11 @@ end
     for algo in algos
         @testset "$algo" begin
             # Extract isosurface using a function
-            points_fn, faces_fn = isosurface(torus_function, algo(), origin=SVector(-2, -2, -2), widths=SVector(4, 4, 4), samples=(81, 81, 81))
+            points_fn, faces_fn = isosurface(torus_function, algo(), -2:2, -2:2, -2:2, samples=(81, 81, 81))
 
             # Extract isosurface using an array
-            torus_array = [torus_function(SVector(x, y, z)) for x in -2:0.05:2, y in -2:0.05:2, z in -2:0.05:2]
-            points_arr, faces_arr = isosurface(torus_array, algo(), origin=SVector(-2, -2, -2), widths=SVector(4, 4, 4))
+            torus_array = [torus_function((x, y, z)) for x in -2:0.05:2, y in -2:0.05:2, z in -2:0.05:2]
+            points_arr, faces_arr = isosurface(torus_array, algo(), -2:2, -2:2, -2:2)
 
             # Test that the vertices and faces are the same for both cases
             @test_broken all(points_fn .≈ points_arr)
@@ -179,17 +130,20 @@ end
     samples = randn(10, 10, 10)
 
     # Extract isosurfaces using MarchingTetrahedra
-    points_mt1, faces_mt1 = isosurface(samples, MarchingTetrahedra(), origin=SVector(0, 0, 0), widths=SVector(1, 1, 1))
-    points_mt2, faces_mt2 = isosurface(Float32.(samples), MarchingTetrahedra(), origin=SVector(0, 0, 0), widths=SVector(1, 1, 1))
-
-    @test all(isapprox.(points_mt1, points_mt2, rtol=1e-6))
-    @test all(faces_mt1 == faces_mt2)
+    points_mt1, faces_mt1 = isosurface(samples, MarchingTetrahedra(), 0:1, 0:1, 0:1)
+    points_mt2, faces_mt2 = isosurface(Float32.(samples), MarchingTetrahedra(), 0:1, 0:1, 0:1)
+
+    @test length(points_mt1) == length(points_mt2)
+    #for i in eachindex(points_mt1)
+    #    @test all(points_mt1[i] .≈ points_mt2[i])
+    #end
+    @test length(faces_mt1) == length(faces_mt2)
+    for i in eachindex(faces_mt1)
+        @test all(faces_mt1[i] .≈ faces_mt2[i])
+    end
 
-    #=
     @testset "forward diff" begin
         # Demonstrate that we can even take derivatives through the meshing algorithm
-        origin = SVector(-1, -1, -1)
-        widths = SVector(2, 2, 2)
         resolution = 0.1
 
         for algo in (MarchingCubes, MarchingTetrahedra)
@@ -203,7 +157,7 @@ end
                         # of some quantity you might want to differentiate in a mesh,
                         # and has the benefit for testing of having a trivial expected
                         # solution.
-                        points, _ = isosurface(norm, algo(isovalue,  reduce_vert=false), samples=(21, 21, 21))
+                        points, _ = isosurface(norm, algo(iso=isovalue), samples=(21, 21, 21))
                         mean(norm.(points))
                     end
 
@@ -215,7 +169,7 @@ end
         end
     end
 
-    =#
+
 
     @testset "type stability" begin
         # verify that we don't lose type stability just by mixing arguments
@@ -223,36 +177,38 @@ end
         data = randn(5, 5, 5)
         iso = 0.2
         eps = 1e-4
-        for algo_ty in (MarchingTetrahedra, MarchingCubes,NaiveSurfaceNets)
-            algo1 = algo_ty(iso, eps)
-            algo2 = algo_ty(Float64(iso), Float16(eps))
-            algo3 = algo_ty(Float32(iso), Float64(eps))
-            @inferred(isosurface(data, algo1, SVector{3, Float16}, SVector{3,UInt32}))
-            @inferred(isosurface(Float32.(data), algo2, SVector{3, Float32}, SVector{3,Int16}))
-            @inferred(isosurface(Float64.(data), algo3, SVector{3, Float64}, SVector{3,UInt}))
-        end
+        algo1 = MarchingTetrahedra(iso=iso, eps=eps)
+        algo2 = MarchingTetrahedra(iso=Float16(iso), eps=Float16(eps))
+        algo3 = MarchingTetrahedra(iso=Float32(iso), eps=Float32(eps))
+        @inferred(isosurface(data, algo1))
+        @inferred(isosurface(Float32.(data), algo2))
+        @inferred(isosurface(Float64.(data), algo3))
+
+        algo1 = MarchingCubes(iso=iso)
+        algo2 = MarchingCubes(iso=Float16(iso))
+        algo3 = MarchingCubes(iso=Float32(iso))
+        @inferred(isosurface(data, algo1))
+        @inferred(isosurface(Float32.(data), algo2))
+        @inferred(isosurface(Float64.(data), algo3))
     end
 
     @testset "Float16" begin
         sdf_torus = [((sqrt(x^2 + y^2) - 0.5)^2 + z^2 - 0.25) for x in -2:0.05:2, y in -2:0.05:2, z in -2:0.05:2]
-        points_nsn, faces_nsn = isosurface(Float16.(sdf_torus), NaiveSurfaceNets(Float16(0.0)), origin=SVector(Float16(-2), Float16(-2), Float16(-2)), widths=SVector(Float16(4), Float16(4), Float16(4)))
-        points_mt, faces_mt = isosurface(Float16.(sdf_torus), MarchingTetrahedra(Float16(0.0)), origin=SVector(Float16(-2), Float16(-2), Float16(-2)), widths=SVector(Float16(4), Float16(4), Float16(4)))
-        points_mc, faces_mc = isosurface(Float16.(sdf_torus), MarchingCubes(Float16(0.0)), origin=SVector(Float16(-2), Float16(-2), Float16(-2)), widths=SVector(Float16(4), Float16(4), Float16(4)))
-
-        @test_broken typeof(points_nsn) == Vector{SVector{3, Float16}}
-        @test typeof(faces_nsn) == Vector{SVector{4, Int}}
-        @test_broken typeof(points_mt) == Vector{SVector{3, Float16}}
-        @test typeof(faces_mt) == Vector{SVector{3, Int}}
-        @test_broken typeof(points_mc) == Vector{SVector{3, Float16}}
-        @test typeof(faces_mc) == Vector{SVector{3, Int}}
+        points_mt, faces_mt = isosurface(Float16.(sdf_torus), MarchingTetrahedra(iso=Float16(0.0), eps=Float16(1e-3)), Float16(-2):Float16(2), Float16(-2):Float16(2), Float16(-2):Float16(2))
+        points_mc, faces_mc = isosurface(Float16.(sdf_torus), MarchingCubes(iso=Float16(0.0)), Float16(-2):Float16(2), Float16(-2):Float16(2), Float16(-2):Float16(2))
+
+        @test typeof(points_mt) == Vector{NTuple{3, Float16}}
+        @test typeof(faces_mt) == Vector{NTuple{3, Int}}
+        @test typeof(points_mc) == Vector{NTuple{3, Float16}}
+        @test typeof(faces_mc) == Vector{NTuple{3, Int}}
     end
     @testset "Integers" begin
         A = rand(Int, 10,10,10)
         for algo in algos
             @testset "$algo" begin
                 p, t = isosurface(A, algo())
-                @test typeof(p) <: Vector{SVector{3,Float64}}
-                @test typeof(t) <: Vector{ SVector{ Meshing.default_face_length(algo()), Int } }
+                @test typeof(p) <: Vector{NTuple{3, Float64}}
+                @test typeof(t) <: Vector{NTuple{3, Int}}
             end
         end
     end