diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 4cef005..693c136 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -3,9 +3,11 @@ on:
   pull_request:
     branches:
       - master
+      - dev
   push:
     branches:
       - master
+      - dev
     tags: '*'
 jobs:
   test:
diff --git a/src/AlphaShapes.jl b/src/AlphaShapes.jl
index f1d489b..615429b 100644
--- a/src/AlphaShapes.jl
+++ b/src/AlphaShapes.jl
@@ -1,7 +1,7 @@
 module AlphaShapes
     import BlackBoxOptim.bboptimize, BlackBoxOptim.best_candidate
     import Distances.pairwise, Distances.Euclidean
-    import LinearAlgebra.det, LinearAlgebra.inv
+    import LinearAlgebra.det, LinearAlgebra.inv, LinearAlgebra.norm
 
     using MiniQhull
 
@@ -123,7 +123,7 @@ julia>:([0    1      1      1
     Wrap MiniQhull.jl's delaunay to get a delaunay triangualation in any
     dimension
     """
-    function GetDelaunayTriangulation(points::AbstractArray{Float64,2})::AbstractArray{Float64,3}
+    function GetDelaunayTriangulation(points::AbstractArray{Float64,2}, indices::Bool=false)
         tess = delaunay(points)
         Triangles = zeros(size(tess,1)-1,size(tess,1),size(tess,2))
         for i in axes(tess, 1)
@@ -131,7 +131,7 @@ julia>:([0    1      1      1
                 view(Triangles, :, i, j) .= view(points, :, tess[i, j])
             end
         end
-        return Triangles
+        return indices ? (Triangles, tess) : Triangles
     end
 
     """
@@ -198,4 +198,7 @@ julia>:([0    1      1      1
         end
         return T[:,:,A]
     end
+
+	include("utils.jl")
+	include("density.jl")
 end # module AlphaShapes
diff --git a/src/density.jl b/src/density.jl
new file mode 100644
index 0000000..68d3338
--- /dev/null
+++ b/src/density.jl
@@ -0,0 +1,79 @@
+"""
+    ContiguousVoronoiCell(point_index::Int,tess,tess_inds)
+
+For a data point (point_index) return the simplices which share this
+point as a common vertex
+"""
+function ContiguousVoronoiCell(point_index::Int,tess,tess_inds)
+    tris = []
+    pos = []
+    for ind in findall(x->x.==point_index,tess_inds)
+        push!(tris,ind[2])
+        push!(pos,ind[1])
+    end
+    return tess[tris,:,:],pos
+end
+
+"""
+    ContiguousVoronoiCellArea(point_index::Int,tess,tess_inds)
+
+For a data point (point_index) return the sum of the areas of the simplices
+which share this point as a common vertex
+"""
+function ContiguousVoronoiCellArea(point_index::Int,tess,tess_inds)
+    tris,inds = ContiguousVoronoiCell(point_index,tess,tess_inds)
+    A = 0.0
+    for i in 1:size(tris,1)
+        A += AlphaShapes.SimplexVolume(tris[i,:,:])
+    end
+    return A
+end
+
+"""
+    DTFELocalDensity(point_index::Int,tess,tess_inds)
+
+Return the local desnity at a point in the Delaunay Tesselation Field
+Estimator sense.
+
+https://arxiv.org/pdf/0708.1441.pdf
+"""
+DTFELocalDensity(point_index::Int,tess,tess_inds)::Float64 = size(tess,2)/ContiguousVoronoiCellArea(point_index,tess,tess_inds)
+
+function LawOfCosine(a::Float64,b::Float64,c::Float64)::Float64
+    return acos(max(-1.0,min(1.0,(a^2 + b^2 - c^2) / (2*a*b+1e-100))))
+end
+
+function TriangleAngles(triangle::Array{Float64,2})::Array{Float64,1}
+    p_a, p_b, p_c = triangle[1,:],triangle[2,:],triangle[3,:]
+    a = norm(p_c .- p_b)
+    b = norm(p_a .- p_c)
+    c = norm(p_a .- p_b)
+
+    A = LawOfCosine(b,c,a)
+    B = LawOfCosine(a,c,b)
+    C = π - (A + B)
+    return [A,B,C]
+end
+
+PointWeights(triangle::Array{Float64,2})::Array{Float64,1} =  TriangleAngles(triangle) ./ 2π
+
+"""
+    WeightedDTFELocalDensity(point_index::Int,tess,tess_inds)
+
+Return the local desnity at a point in the Weighted Delaunay Tesselation Field
+Estimator sense. This additionally normalises by the 'contribution' of individual
+points (portion of 2pi rad. in triangle)
+
+ https://royalsocietypublishing.org/doi/10.1098/rsif.2021.0114
+"""
+
+function WeightedDTFELocalDensity(point_index::Int,tess,tess_inds)::Float64
+    tris,pos = ContiguousVoronoiCell(point_index,tess,tess_inds)
+    A = zeros(size(tris,1))
+    w = zeros(size(tris,1))
+    for i in 1:size(tris,1)
+        w[i] = PointWeights(tris[i,:,:])[pos[i]]
+        A[i] = AlphaShapes.SimplexVolume(tris[i,:,:])
+    end
+    return sum(w./A)/2.0*π
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 1d10012..c593386 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,7 +2,8 @@ using Test
 my_tests = [
   "test_circum_sphere.jl",
   "test_simplex_volume.jl",
-  "test_alpha_shape.jl"
+  "test_alpha_shape.jl",
+  "test_density.jl"
 ]
 
 println("Running tests:")
diff --git a/test/test_density.jl b/test/test_density.jl
new file mode 100644
index 0000000..cf84b5f
--- /dev/null
+++ b/test/test_density.jl
@@ -0,0 +1,25 @@
+using AlphaShapes
+using Test
+
+@testset "Test Local Density" begin
+    # four corners
+    X = [0.0 0.0 1.0 1.0; 0.0 1.0 1.0 0.0]
+    tess, inds = AlphaShapes.GetDelaunayTriangulation(X,true)
+    density = [AlphaShapes.WeightedDTFELocalDensity(i,tess,inds) for i in 1:size(X,2)]
+
+    @test density == [
+        0.39269908169872425,
+        0.39269908169872403,
+        0.39269908169872425,
+        0.39269908169872403
+    ]
+
+    density = [AlphaShapes.DTFELocalDensity(i,tess,inds) for i in 1:size(X,2)]
+
+    @test density == [
+        3.0,
+        1.5,
+        3.0,
+        1.5       
+    ]
+end
\ No newline at end of file