diff --git a/src/Meshes.jl b/src/Meshes.jl
index 2bef7cd51..723373596 100644
--- a/src/Meshes.jl
+++ b/src/Meshes.jl
@@ -519,6 +519,7 @@ export
   JarvisMarch,
   hull,
   convexhull,
+  concavehull,
 
   # sampling
   SamplingMethod,
diff --git a/src/hulls.jl b/src/hulls.jl
index dc11291d8..a664c98b4 100644
--- a/src/hulls.jl
+++ b/src/hulls.jl
@@ -34,6 +34,13 @@ Convex hull of `object`.
 """
 function convexhull end
 
+"""
+    concavehull(object)
+
+Concave hull of `object`.
+"""
+# function concavehull end
+
 # ----------
 # FALLBACKS
 # ----------
@@ -46,6 +53,14 @@ convexhull(m::Multi) = _gconvexhull(parent(m))
 
 convexhull(geoms) = _gconvexhull(geoms)
 
+concavehull(p::Polytope) = _pconcavehull(eachvertex(p))
+
+concavehull(p::Primitive) = concavehull(boundary(p))
+
+concavehull(m::Multi) = _gconcavehull(parent(m))
+
+concavehull(geoms) = _gconcavehull(geoms)
+
 # ----------------
 # SPECIALIZATIONS
 # ----------------
@@ -64,6 +79,20 @@ convexhull(g::Grid) = Box(extrema(g)...)
 
 convexhull(m::Mesh) = _pconvexhull(eachvertex(m))
 
+concavehull(p::Point) = p
+
+concavehull(b::Box) = b
+
+concavehull(b::Ball) = b
+
+concavehull(s::Sphere) = Ball(center(s), radius(s))
+
+concavehull(t::Triangle) = t
+
+concavehull(g::Grid) = Box(extrema(g)...)
+
+concavehull(m::Mesh) = _pconcavehull(eachvertex(m))
+
 # ----------------
 # IMPLEMENTATIONS
 # ----------------
@@ -71,3 +100,7 @@ convexhull(m::Mesh) = _pconvexhull(eachvertex(m))
 _gconvexhull(geoms) = _pconvexhull(p for g in geoms for p in boundarypoints(g))
 
 _pconvexhull(points) = hull(points, GrahamScan())
+
+_gconcavehull(geoms) = _pconcavehull(p for g in geoms for p in boundarypoints(g))
+
+_pconcavehull(points) = hull(points, JarvisMarch(3)) # default k = 3 for concave hulls
diff --git a/src/hulls/jarvis.jl b/src/hulls/jarvis.jl
index 0bdcc4397..ca45f1e5e 100644
--- a/src/hulls/jarvis.jl
+++ b/src/hulls/jarvis.jl
@@ -4,11 +4,15 @@
 
 """
     JarvisMarch()
+    JarvisMarch(k)
 
 Compute the convex hull of a set of points or geometries using the
 Jarvis's march algorithm. See [https://en.wikipedia.org/wiki/Gift_wrapping_algorithm]
 (https://en.wikipedia.org/wiki/Gift_wrapping_algorithm).
 
+If `k` is provided, the algorithm will attempt to compute a concave hull using [k,n)
+nearest neighbors.
+
 The algorithm has complexity `O(n*h)` where `n` is the number of points
 and `h` is the number of points in the hull.
 
@@ -16,10 +20,15 @@ and `h` is the number of points in the hull.
 
 * Jarvis 1973. [On the identification of the convex hull of a finite set of
   points in the plane](https://www.sciencedirect.com/science/article/abs/pii/0020019073900203)
+* Moreira, A. & Santos, M. Y. 2007. "concave hull: a k-nearest neighbours approach for the computation of the region occupied by a set of points". In: *Proceedings of the Second International Conference on Computer Graphics Theory and Applications*. SciTePress, pp. 61–68.
 """
-struct JarvisMarch <: HullMethod end
+struct JarvisMarch{K} <: HullMethod
+  k::K
+end
+JarvisMarch() = JarvisMarch{Nothing}(nothing)
+JarvisMarch(k::I) where {I<:Integer} = JarvisMarch{I}(k)
 
-function hull(points, ::JarvisMarch)
+function hull(points, method::JarvisMarch)
   pₒ = first(points)
   ℒ = lentype(pₒ)
 
@@ -34,16 +43,29 @@ function hull(points, ::JarvisMarch)
   # corner cases
   n == 1 && return p[1]
   n == 2 && return Segment(p[1], p[2])
+  n == 3 && return PolyArea(p)
 
   # find bottom-left point
   i = argmin(p)
+  # find next point with smallest angle
+  O = p[i]
+  A = O + Vec(zero(ℒ), -oneunit(ℒ))
+
+  # perform primary Jarvis march loop
+  ℐ = _jarvisloop(method, points, p, n, i, O, A)
+  # if no hull is found, return convex hull as fallback
+  isnothing(ℐ) && return convexhull(p)
+
+  # return polygonal area
+  PolyArea(p[ℐ[begin:(end - 1)]])
+end
 
+# convex hull
+function _jarvisloop(::JarvisMarch{Nothing}, points, p, n, i, O, A)
   # candidates for next point
   𝒞 = [1:(i - 1); (i + 1):n]
 
   # find next point with smallest angle
-  O = p[i]
-  A = O + Vec(zero(ℒ), -oneunit(ℒ))
   j = argmin(l -> ∠(A, O, p[l]), 𝒞)
 
   # initialize ring of indices
@@ -67,6 +89,92 @@ function hull(points, ::JarvisMarch)
     push!(ℐ, j)
   end
 
-  # return polygonal area
-  PolyArea(p[ℐ[begin:(end - 1)]])
+  # return indices of hull vertices
+  ℐ
+end
+
+# concave hull
+function _jarvisloop(method::JarvisMarch{I}, points, p, n, i, O, A) where {I<:Integer}
+  m = I(n - 2)
+  k = min(max(method.k, 3), m)
+  assertion(ispositive(k), "k must be a positive integer")
+  k > m && return _jarvisloop(JarvisMarch{Nothing}(nothing), points, p, n, i, O, A) # fallback to convex hull
+
+  # initial state for retries
+  i₀, O₀, A₀ = i, O, A
+
+  # try increasing k until valid hull found
+  for ki in k:m
+    searcher = KNearestSearch(p, ki)
+    mask = trues(n)
+
+    # apply initial point information
+    i = i₀
+    O = O₀
+    A = A₀
+    mask[i] = false
+
+    # find next point with smallest angle among k nearest neighbors
+    𝒩 = search(O, searcher; mask=mask)
+    isempty(𝒩) && (k += 1; continue)
+
+    j = argmin(l -> ∠(A, O, p[l]), 𝒩)
+    ℐ = [i, j]
+
+    # no longer searchable
+    mask[j] = false
+    step = 2
+    failed = false
+
+    while first(ℐ) != last(ℐ)
+      # reinsert first point after 5 steps. Otherwise the searcher may double back and fail to find a valid hull.
+      step == 5 && (mask[ℐ[begin]] = true)
+
+      # direction of current segment
+      v = p[j] - p[i]
+      i = j
+      O = p[i]
+      A = O + v
+
+      # order neighbors by angle
+      search!(𝒩, O, searcher; mask)
+      isempty(𝒩) && break
+      sort!(𝒩, by=l -> ∠(A, O, p[l]))
+
+      found = false
+
+      for nᵢ in 𝒩
+        cpoint = p[nᵢ]
+        last = cpoint == p[ℐ[begin]] ? 1 : 0
+
+        its = false
+        indⱼ = 2
+        while !its && indⱼ < length(ℐ) - last
+          its = intersects(Segment(p[ℐ[step]], cpoint), Segment(p[ℐ[step - indⱼ + 1]], p[ℐ[step - indⱼ]]))
+          indⱼ += 1
+        end
+
+        if !its
+          j = nᵢ
+          found = true
+          break
+        end
+      end
+
+      # if no candidate found, break and try again with larger k
+      !found && (failed = true; break)
+
+      mask[j] = false
+      step += 1
+      push!(ℐ, j)
+    end
+
+    # check if hull is valid
+    if !failed
+      poly = PolyArea(p[ℐ[begin:(end - 1)]])
+      all(points .∈ poly) && return ℐ
+    end
+  end
+
+  nothing
 end
diff --git a/test/hulls.jl b/test/hulls.jl
index 1717a4eee..d0135465b 100644
--- a/test/hulls.jl
+++ b/test/hulls.jl
@@ -1,5 +1,5 @@
 @testitem "Hulls" setup = [Setup] begin
-  for method in [GrahamScan(), JarvisMarch()]
+  for method in [GrahamScan(), JarvisMarch(), JarvisMarch(3)]
     # basic test
     pts = [cart(rand(T), rand(T)) for _ in 1:10]
     chul = hull(pts, method)
@@ -20,21 +20,25 @@
     @test chul == Segment(cart(0, 1), cart(1, 0))
     pts = cart.([(1, 0), (0, 0), (0, 1)])
     chul = hull(pts, method)
-    @test vertices(chul) == cart.([(0, 0), (1, 0), (0, 1)])
+    @test Set(vertices(chul)) == Set(cart.([(0, 0), (1, 0), (0, 1)]))
 
     # original point set is already in hull
     pts = cart.([(0, 0), (1, 0), (1, 1), (0, 1), (0.5, -1)])
     chul = hull(pts, method)
     verts = vertices(chul)
-    @test verts == cart.([(0, 0), (0.5, -1), (1, 0), (1, 1), (0, 1)])
+    @test Set(verts) == Set(cart.([(0, 0), (0.5, -1), (1, 0), (1, 1), (0, 1)]))
 
-    # random points in interior do not affect result
+    # all points should be in hull, even if random 
     p1 = cart.([(0, 0), (1, 0), (1, 1), (0, 1), (0.5, -1)])
     p2 = cart.([0.5 .* (rand(), rand()) .+ 0.5 for _ in 1:10])
     pts = [p1; p2]
     chul = hull(pts, method)
     verts = vertices(chul)
-    @test verts == cart.([(0, 0), (0.5, -1), (1, 0), (1, 1), (0, 1)])
+    @test all(p1 .∈ Ref(chul))
+    @test all(p2 .∈ Ref(chul))
+    if !(method isa JarvisMarch && !isnothing(method.k)) # convex hull should be unaffected by interior points
+      @test verts == cart.([(0, 0), (0.5, -1), (1, 0), (1, 1), (0, 1)])
+    end
 
     pts =
       cart.([
@@ -66,12 +70,20 @@
         (0, 6)
       ])
     chul = hull(pts, method)
-    @test nvertices(chul) < length(pts)
+    if method isa JarvisMarch && !isnothing(method.k)
+      @test Set(vertices(chul)) == Set(pts)
+    else
+      @test nvertices(chul) < length(pts)
+    end
 
     poly = readpoly(T, joinpath(datadir, "hull.line"))
     pts = vertices(poly)
     chul = hull(pts, method)
-    @test nvertices(chul) < length(pts)
+    if method isa JarvisMarch && !isnothing(method.k)
+      @test Set(vertices(chul)) == Set(pts)
+    else
+      @test nvertices(chul) < length(pts)
+    end
 
     if method == GrahamScan()
       # simplifying rectangular hull / triangular
@@ -171,3 +183,26 @@ end
   @test cart(-0.8, -0.8) ∈ h
   @test cart(0.2, 0.2) ∈ h
 end
+
+@testitem "Concave hulls" setup = [Setup] begin
+  @test concavehull(cart(0, 0)) == cart(0, 0)
+
+  @test concavehull(Box(cart(0, 0), cart(1, 1))) == Box(cart(0, 0), cart(1, 1))
+
+  @test concavehull(Ball(cart(0, 0), T(1))) == Ball(cart(0, 0), T(1))
+  @test concavehull(Ball(cart(1, 1), T(1))) == Ball(cart(1, 1), T(1))
+
+  @test concavehull(Sphere(cart(0, 0), T(1))) == Ball(cart(0, 0), T(1))
+  @test concavehull(Sphere(cart(1, 1), T(1))) == Ball(cart(1, 1), T(1))
+
+  b1 = Box(cart(0, 0), cart(1, 1))
+  b2 = Box(cart(-1, -1), cart(0.5, 0.5))
+  @test Set(vertices(concavehull(Multi([b1, b2])))) ==
+        Set(cart.([(-1, -1), (0.5, -1), (1, 0), (1, 1), (0, 1), (-1, 0.5)]))
+
+  b1 = Ball(cart(0, 0), T(1))
+  b2 = Box(cart(-1, -1), cart(0.5, 0.5))
+  h = concavehull(Multi([b1, b2]))
+  @test h isa PolyArea
+  @test all(Set([boundarypoints(b1); boundarypoints(b2)]) .∈ Ref(h))
+end