Fast concave hull implementation in Python.

ConcaveHull.py

'''
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__':

how to return the boundary points?

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.

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

This code produces CONVEX hull:

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

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.

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.

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.

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

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.

@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?

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

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