Skip to content

Instantly share code, notes, and snippets.

@AndreLester
Last active January 16, 2024 18:15
  • Star 17 You must be signed in to star a gist
  • Fork 2 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save AndreLester/589ea1eddd3a28d00f3d7e47bd9f28fb to your computer and use it in GitHub Desktop.
Fast concave hull implementation in Python.
'''
Copyright (C) 2018 Andre Lester Kruger
ConcaveHull.py is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 2 of the License, or
(at your option) any later version.
ConcaveHull.py is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ConcaveHull.py. If not, see <http://www.gnu.org/licenses/>.
'''
import bisect
from collections import OrderedDict
import math
#import numpy as np
import matplotlib.tri as tri
from shapely.geometry import LineString
from shapely.geometry import Polygon
from shapely.ops import linemerge
class ConcaveHull:
def __init__(self):
self.triangles = {}
self.crs = {}
def loadpoints(self, points):
#self.points = np.array(points)
self.points = points
def edge(self, key, triangle):
'''Calculate the length of the triangle's outside edge
and returns the [length, key]'''
pos = triangle[1].index(-1)
if pos==0:
x1, y1 = self.points[triangle[0][0]]
x2, y2 = self.points[triangle[0][1]]
elif pos==1:
x1, y1 = self.points[triangle[0][1]]
x2, y2 = self.points[triangle[0][2]]
elif pos==2:
x1, y1 = self.points[triangle[0][0]]
x2, y2 = self.points[triangle[0][2]]
length = ((x1-x2)**2+(y1-y2)**2)**0.5
rec = [length, key]
return rec
def triangulate(self):
if len(self.points) < 2:
raise Exception('CountError: You need at least 3 points to Triangulate')
temp = list(zip(*self.points))
x, y = list(temp[0]), list(temp[1])
del(temp)
triang = tri.Triangulation(x, y)
self.triangles = {}
for i, triangle in enumerate(triang.triangles):
self.triangles[i] = [list(triangle), list(triang.neighbors[i])]
def calculatehull(self, tol=50):
self.tol = tol
if len(self.triangles) == 0:
self.triangulate()
# All triangles with one boundary longer than the tolerance (self.tol)
# is added to a sorted deletion list.
# The list is kept sorted from according to the boundary edge's length
# using bisect
deletion = []
self.boundary_vertices = set()
for i, triangle in self.triangles.items():
if -1 in triangle[1]:
for pos, neigh in enumerate(triangle[1]):
if neigh == -1:
if pos == 0:
self.boundary_vertices.add(triangle[0][0])
self.boundary_vertices.add(triangle[0][1])
elif pos == 1:
self.boundary_vertices.add(triangle[0][1])
self.boundary_vertices.add(triangle[0][2])
elif pos == 2:
self.boundary_vertices.add(triangle[0][0])
self.boundary_vertices.add(triangle[0][2])
if -1 in triangle[1] and triangle[1].count(-1) == 1:
rec = self.edge(i, triangle)
if rec[0] > self.tol and triangle[1].count(-1) == 1:
bisect.insort(deletion, rec)
while len(deletion) != 0:
# The triangles with the longest boundary edges will be
# deleted first
item = deletion.pop()
ref = item[1]
flag = 0
# Triangle will not be deleted if it already has two boundary edges
if self.triangles[ref][1].count(-1) > 1:
continue
# Triangle will not be deleted if the inside node which is not
# on this triangle's boundary is already on the boundary of
# another triangle
adjust = {0: 2, 1: 0, 2: 1}
for i, neigh in enumerate(self.triangles[ref][1]):
j = adjust[i]
if neigh == -1 and self.triangles[ref][0][j] in self.boundary_vertices:
flag = 1
break
if flag == 1:
continue
for i, neigh in enumerate(self.triangles[ref][1]):
if neigh == -1:
continue
pos = self.triangles[neigh][1].index(ref)
self.triangles[neigh][1][pos] = -1
rec = self.edge(neigh, self.triangles[neigh])
if rec[0] > self.tol and self.triangles[rec[1]][1].count(-1) == 1:
bisect.insort(deletion, rec)
for pt in self.triangles[ref][0]:
self.boundary_vertices.add(pt)
del self.triangles[ref]
self.polygon()
def polygon(self):
edgelines = []
for i, triangle in self.triangles.items():
if -1 in triangle[1]:
for pos, value in enumerate(triangle[1]):
if value == -1:
if pos==0:
x1, y1 = self.points[triangle[0][0]]
x2, y2 = self.points[triangle[0][1]]
elif pos==1:
x1, y1 = self.points[triangle[0][1]]
x2, y2 = self.points[triangle[0][2]]
elif pos==2:
x1, y1 = self.points[triangle[0][0]]
x2, y2 = self.points[triangle[0][2]]
line = LineString([(x1, y1), (x2, y2)])
edgelines.append(line)
bound = linemerge(edgelines)
self.boundary = Polygon(bound.coords)
#if __name__ == '__main__':
@RobSimpsonHove
Copy link

Thanks, I thought as much. I've ended up using alphashape.

@TinKurbatoff
Copy link

TinKurbatoff commented Apr 25, 2021

I am using alpha shape, but it still not as well as I want to:
image

——————
It is a "manual" code, and it used the Delaunay function.
The code template I took from here:
https://deeplearning.lipingyang.org/wp-content/uploads/2019/07/Drawing-Boundaries-In-Python.pdf
Not ideal, but it works better.
—————
I still try to build a concave hull algorithm implementation. Will keep you posted if succeed.

@arjun-bala-krishnan
Copy link

I am still trying to find a similar implementation in 3D using python. If anybody manages to find it (similar to MATLAB boundary function for 3d) it would be really helpful.

@jatinnimawat
Copy link

Is this Concave or Convex Hulls - I am getting only Convex hulls from this.

I'm also getting Convex hulls instead of the intended Concave ones. This code here provided better results, at least in my case, although they are not yet what I need. Hope it helps you.

https://gist.github.com/dwyerk/10561690

Same here. I am also getting convex hulls instead of concave.

@cvr
Copy link

cvr commented Feb 8, 2023

To get a concave hull you need to change the tolerance by changing argument tol:

ch = ConcaveHull()
ch.loadpoints(pts)
ch.calculatehull(tol=0.5)

Is this Concave or Convex Hulls - I am getting only Convex hulls from this.

I'm also getting Convex hulls instead of the intended Concave ones. This code here provided better results, at least in my case, although they are not yet what I need. Hope it helps you.
https://gist.github.com/dwyerk/10561690

Same here. I am also getting convex hulls instead of concave.

@AndreLester
Copy link
Author

AndreLester commented Feb 8, 2023

@cvr The "tol" variable presents the length of the longest edge you wish to have in the respective units. Usually in feet or metres. It is not a factor between zero and one as other concave hull algorithms. I prefer it like this.
The problems above where the people got convex hulls was probably due to the fact that their units was in degrees.

Note: Apologies. I was unaware that my "gist" blew up like this.

@cvr
Copy link

cvr commented Feb 8, 2023

Thanks for the clarification. Matter of fact I was also using degrees in the test case where I'm trying your code.

It actually makes a lot of sense now. You could have it in degrees, it is just a question of changing the tol value to a meaningful value of distance, if you would convert it to degrees (e.g. divide tol in km it by 111.045 to have a rough measure in degrees).

As this will depend on some triangulation of a set of points, eventually tol could be made as a fraction of the total width (or height) span of the collection of points.

@Joshua-Mursic
Copy link

@AndreLester, Thanks a lot for the code! I am trying to implement a offset value into this function, do you have any idea how this would work?
example

Like in the image above the green is original shape and the blue is the offset shape.

@AndreLester
Copy link
Author

AndreLester commented Jun 2, 2023

@Joshua-Mursic I would use the buffer function from shapely on the resultant polygon. Play around with the options to get the desired result.

@AndreLester
Copy link
Author

I have updated, modified and documented this script. It can be found at https://github.com/civildot/cdBoundary and installed with:

pip install cdBoundary

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment