''' | |
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__': |
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
Much faster then using alphashape!
how can i plot the surface of the polygon?
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
Used your code as part of my project. Thought it would interesting for you:
https://externalflow.et.aau.dk/image-processing-techniques/
Is this Concave or Convex Hulls - I am getting only Convex hulls from this.
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.
Thanks, I thought as much. I've ended up using alphashape.
I am using alpha shape, but it still not as well as I want to:
——————
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.
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.
how to return the boundary points?