diff --git a/Project.toml b/Project.toml
index afdc390..fd1851a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -2,7 +2,11 @@ name = "Meshing"
 uuid = "e6723b4c-ebff-59f1-b4b7-d97aa5274f73"
 
 [deps]
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GeometryTypes = "4d00f742-c7ba-57c2-abde-4428a4b178cb"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [extras]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
diff --git a/README.md b/README.md
index 020e439..40f00d1 100644
--- a/README.md
+++ b/README.md
@@ -9,6 +9,7 @@ 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/)
+* [Dual Contours](https://en.wikipedia.org/wiki/Isosurface#Dual_Contouring)
 
 ## Interface
 
@@ -44,6 +45,7 @@ Three meshing algorithms exist:
 * `MarchingCubes()`
 * `MarchingTetrahedra()`
 * `NaiveSurfaceNets()`
+* `DualContours()`
 
 Each takes an optional `iso` and `eps` parameter, e.g. `MarchingCubes(0.0,1e-6)`.
 
diff --git a/src/Meshing.jl b/src/Meshing.jl
index a8218a6..d0bfae6 100644
--- a/src/Meshing.jl
+++ b/src/Meshing.jl
@@ -1,16 +1,22 @@
 module Meshing
 
 using GeometryTypes
+using Roots
+using ForwardDiff
+using LinearAlgebra
+using StaticArrays
 
 abstract type AbstractMeshingAlgorithm end
 
 include("marching_tetrahedra.jl")
 include("marching_cubes.jl")
 include("surface_nets.jl")
+include("dual_contours.jl")
 
 export marching_cubes,
        MarchingCubes,
        MarchingTetrahedra,
-       NaiveSurfaceNets
+       NaiveSurfaceNets,
+       DualContours
 
 end # module
diff --git a/src/dual_contours.jl b/src/dual_contours.jl
new file mode 100644
index 0000000..88b42da
--- /dev/null
+++ b/src/dual_contours.jl
@@ -0,0 +1,128 @@
+#Cardinal directions
+const dirs = ( Point(1,0,0), Point(0,1,0), Point(0,0,1) )
+
+#Vertices of cube
+const cube_verts = (Point(0, 0, 0), Point(0, 0, 1), Point(0, 1, 0),
+                    Point(0, 1, 1), Point(1, 0, 0), Point(1, 0, 1),
+                    Point(1, 1, 0), Point(1, 1, 1))
+
+
+#Edges of cube
+const cube_edges_dc = ((0, 1), (0, 2), (0, 1), (0, 4), (0, 2), (0, 4), (2, 3), (1, 3),
+              (4, 5), (1, 5), (4, 6), (2, 6), (4, 5), (4, 6), (2, 3), (2, 6),
+              (1, 3), (1, 5), (6, 7), (5, 7), (6, 7), (3, 7), (5, 7), (3, 7))
+
+
+#Use non-linear root finding to compute intersection point
+function estimate_hermite(f, v0, v1)
+    l(t) = f((1.0-t)*v0 + t*v1)
+    dl(t) = ForwardDiff.derivative(l,t)
+    t0 = find_zero((l,dl),(0, 1))
+    x0 = (1.0-t0)*v0 + t0*v1
+    return (x0, ForwardDiff.gradient(f,x0))
+end
+
+
+function dual_contours(f::Function, 
+                       bounds::HyperRectangle,
+                       samples::NTuple{3,Int}=(128,128,128),
+                       iso=0.0,
+                       MT::Type{M}=SimpleMesh{Point{3,Float64},Face{3,Int}},
+                       eps=0.00001) where {ST,FT,M<:AbstractMesh}
+
+    orig = origin(bounds)
+    width = widths(bounds)
+    scale = width ./ Point(samples)
+    #Compute vertices
+    dc_verts = []
+    vindex   = Dict()
+    for x in 0:samples[1], y in 0:samples[2], z in 0:samples[3]
+        idx = Point(x,y,z)
+        o = Point(x,y,z) .* scale + orig
+
+        #Get signs for cube
+        cube_signs = [ f(o+(v.*scale))>0 for v in cube_verts ]
+
+        if all(cube_signs) || !any(cube_signs)
+            continue
+        end
+
+        #Estimate hermite data
+        h_data = [ estimate_hermite(f, o+cube_verts[e[1]+1].*scale, o+cube_verts[e[2]+1].*scale)
+            for e in cube_edges_dc if cube_signs[e[1]+1] != cube_signs[e[2]+1] ]
+
+        #Solve qef to get vertex
+        A = Array{Float64,2}(undef,length(h_data),3)
+        for i in eachindex(h_data)
+            A[i,:] = h_data[i][2]
+        end
+        b = [ dot(d[1],d[2]) for d in h_data ]
+
+        #compute the vertex using pseudo inverse
+        v = pinv(A)*b
+
+        #Throw out failed solutions
+        if norm(v-o) > 2
+            continue
+        end
+
+        #Emit one vertex per every cube that crosses
+        push!(vindex, idx => length(dc_verts))
+        push!(dc_verts, (v, ForwardDiff.gradient(f,v)))
+    end
+
+    #Construct faces
+    dc_faces = Face[]
+    for x in 0:samples[1], y in 0:samples[2], z in 0:samples[3]
+
+        idx = Point(x,y,z)
+        if !haskey(vindex,idx)
+            continue
+        end
+
+        #Emit one face per each edge that crosses
+        o = Point(x,y,z) .* scale + orig
+        for i in (1,2,3)
+            for j in 1:i
+                if haskey(vindex,idx+dirs[i]) && haskey(vindex,idx + dirs[j]) && haskey(vindex,idx + dirs[i] + dirs[j])
+                    # determine orientation of the face from the true normal
+                    v1, tn1 = dc_verts[vindex[idx]+1]
+                    v2, tn2 = dc_verts[vindex[idx+dirs[i]]+1]
+                    v3, tn3 = dc_verts[vindex[idx+dirs[j]]+1]
+
+                    e1 = v1-v2
+                    e2 = v1-v3
+                    c = cross(e1,e2)
+                    if dot(c, tn1) > 0
+                        push!(dc_faces, Face{3,Int}(vindex[idx]+1, vindex[idx+dirs[i]]+1, vindex[idx+dirs[j]]+1) )
+                        push!(dc_faces, Face{3,Int}(vindex[idx+dirs[i]+dirs[j]]+1, vindex[idx+dirs[j]]+1, vindex[idx+dirs[i]]+1) )
+                    else
+                        push!(dc_faces, Face{3,Int}(vindex[idx]+1, vindex[idx+dirs[j]]+1, vindex[idx+dirs[i]]+1) )
+                        push!(dc_faces, Face{3,Int}(vindex[idx+dirs[i]+dirs[j]]+1, vindex[idx+dirs[i]]+1, vindex[idx+dirs[j]]+1) )
+                    end
+                end
+            end
+        end
+
+    end
+    return MT([Point(v[1]...) for v in dc_verts], dc_faces)
+end
+
+struct DualContours{T} <: AbstractMeshingAlgorithm
+    iso::T
+    eps::T
+end
+
+DualContours(iso::T1=0.0, eps::T2=1e-3) where {T1, T2} = DualContours{promote_type(T1, T2)}(iso, eps)
+
+function (::Type{MT})(df::SignedDistanceField, method::DualContours)::MT where {MT <: AbstractMesh}
+    dual_contours(df, method.iso, MT, method.eps)
+end
+
+function (::Type{MT})(f::Function, h::HyperRectangle, size::NTuple{3,Int}, method::DualContours)::MT where {MT <: AbstractMesh}
+    dual_contours(f, h, size, method.iso, MT, method.eps)
+end
+
+function (::Type{MT})(f::Function, h::HyperRectangle, method::DualContours; size::NTuple{3,Int}=(128,128,128))::MT where {MT <: AbstractMesh}
+    dual_contours(f, h, size, method.iso, MT, method.eps)
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index e36573d..e359e8a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -141,6 +141,15 @@ using LinearAlgebra: dot, norm
         @test minimum(vertices(mesh)) ≈ [-0.5, -0.5, -0.5] atol=0.2
     end
 
+    @testset "Dual Contours" begin
+
+        sphere(v) = sqrt(sum(dot(v,v))) - 1 # sphere
+
+        m = SimpleMesh(sphere, HyperRectangle(Vec(-1,-1,-1.),Vec(2,2,2.)), (50,50,50), DualContours())
+        @test length(vertices(m)) == 11754
+        @test length(faces(m)) == 51186
+    end
+
     @testset "AbstractMeshingAlgorithm interface" begin
         f = x -> norm(x) - 0.5
         bounds = HyperRectangle(Vec(-1, -1, -1), Vec(2, 2, 2))