Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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__':
@kinreehou

This comment has been minimized.

Copy link

@kinreehou kinreehou commented Feb 11, 2020

how to return the boundary points?

@jafekb

This comment has been minimized.

Copy link

@jafekb jafekb commented Mar 5, 2020

import numpy as np
from ConcaveHull import ConcaveHull

ch = ConcaveHull()
pts = np.random.uniform(size=(100, 2))
ch.loadpoints(pts)
ch.calculatehull()

boundary_points = np.vstack(ch.boundary.exterior.coords.xy).T
# boundary_points is a subset of pts corresponding to the concave hull

@paulaceccon

This comment has been minimized.

Copy link

@paulaceccon paulaceccon commented Mar 9, 2020

Much faster then using alphashape!

@shiernee

This comment has been minimized.

Copy link

@shiernee shiernee commented Apr 10, 2020

how can i plot the surface of the polygon?

@planzaprecisionag

This comment has been minimized.

Copy link

@planzaprecisionag planzaprecisionag commented Apr 17, 2020

how can i plot the surface of the polygon?

Have you tried accessing the class's polygon object? So something like ploygonToPlot = ch.boundary, then just plot it however you normally plot polygons.

- had typed ch.polygon instead of boundary to reference the polygon object

@jakobhaervig

This comment has been minimized.

Copy link

@jakobhaervig jakobhaervig commented Feb 2, 2021

Used your code as part of my project. Thought it would interesting for you:
https://externalflow.et.aau.dk/image-processing-techniques/

@RobSimpsonHove

This comment has been minimized.

Copy link

@RobSimpsonHove RobSimpsonHove commented Mar 30, 2021

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

@Fajikuma

This comment has been minimized.

Copy link

@Fajikuma Fajikuma commented Apr 21, 2021

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

@TinKurbatoff

This comment has been minimized.

Copy link

@TinKurbatoff TinKurbatoff commented Apr 24, 2021

This code produces CONVEX hull:
image

That is good and clean code, but the naming is incorrect.

@RobSimpsonHove

This comment has been minimized.

Copy link

@RobSimpsonHove RobSimpsonHove commented Apr 25, 2021

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

@TinKurbatoff

This comment has been minimized.

Copy link

@TinKurbatoff 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.

@arjunbkool

This comment has been minimized.

Copy link

@arjunbkool arjunbkool commented Sep 8, 2021

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.

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