Does anyone have experience generating a Concave Hull (Alpha Shape) from a list of points with IronPython? Would greatly appreciate any feedback. Thank you in advance!
1 Like
A quick search brought up this…
A handy couple of nodes being used (which you can dig around in) and some heavyweight package authors 
unfortunately very few people are doing what you need…
Hope that helps,
Mark
1 Like
There’s a concave hull algorithm here (in python) which looks easy to implement.
And there’s some links by @Michael_Kirschner2 here which might offer some alternatives:
Thomas,
I’ve found a couple different concave hull algorithms in python; however, all of them use packages (numpy, scipy, matplotlib, shapely) which are not easily or if at all compatible with Iron Python.
So, I was hoping maybe someone had created a concave hull algorithm with vanilla Python 2.7 modules which could be used in Dynamo.
Mark,
Thank you for the reference. I have reviewed this thread and tried convex hulls for my application but did not find any solution.
Maybe I could reach out to the authors as you suggested. I’m hoping maybe someone has created a concave hull algorithm with vanilla Python 2.7 modules which could be used in Dynamo, but from my research so far there is very few resources regarding this subject.
@ericabbott, have you already go through with Lunchbox ML LunchBoxML for Dynamo - PROVING GROUND
by @Nathan_Miller . Right now, that’s the furthest you can get.
1 Like
Hey,
The link here that MK gave you only uses Math?
There is a C# from string node, but I’m not skilled enough to use it. Or it’s possible you could persuade a package manager to add it…
I had a quick go at translating it to Python, perhaps someone with a maths degree can have a look at it?! But I know @Daniel_Woodcock1 is very busy…
# Enable Python support and load DesignScript library
import clr
clr.AddReference('ProtoGeometry')
from Autodesk.DesignScript.Geometry import *
import math
points = UnwrapElement(IN[0])
#0. error checking, init
if (points == None or points.Count < 2) :
raise Exception('AlphaShape needs at least 2 points')
BorderEdges = []
alpha_2 = alpha * alpha
#1. run through all pairs of points
i = 0
while i < points.Count:
j = i + 1
while j < points.Count:
if (points[i] == points[j]):
raise Exception('AlphaShape needs pairwise distinct points')
dist = (points[i], points[j])
if (dist > 2 * alpha) : continue #circle fits between points ==> p_i, p_j can't be alpha-exposed
x1 = points[i].X
x2 = points[j].X
y1 = points[i].Y
y2 = points[j].Y
mid = PointByCoordinates((x1 + x2) / 2, (y1 + y2) / 2)
#find two circles that contain p_i and p_j; note that center1 == center2 if dist == 2*alpha
center1 = PointByCoordinates((mid.X + math.Sqrt(alpha_2 - (dist / 2) * (dist / 2)) * (y1 - y2) / dist), (mid.Y + math.Sqrt(alpha_2 - (dist / 2) * (dist / 2)) * (x2 - x1) / dist))
center2 = PointByCoordinates((mid.X - math.Sqrt(alpha_2 - (dist / 2) * (dist / 2)) * (y1 - y2) / dist), (mid.Y - math.Sqrt(alpha_2 - (dist / 2) * (dist / 2)) * (x2 - x1) / dist))
#check if one of the circles is alpha-exposed, i.e. no other point lies in it
c1_empty = true
c2_empty = true
k = 0
while k < points.Count & (c1_empty or c2_empty):
if (points[k] == points[i] or points[k] == points[j]):
continue
if ((center1.X - points[k].X) * (center1.X - points[k].X) + (center1.Y - points[k].Y) * (center1.Y - points[k].Y) < alpha_2):
c1_empty = false
if ((center2.X - points[k].X) * (center2.X - points[k].X) + (center2.Y - points[k].Y) * (center2.Y - points[k].Y) < alpha_2):
c2_empty = false
if (c1_empty or c2_empty):
#yup!
BorderEdges.append(Line.ByStartEndPoint(points[i], points[j]))
k += 1
j += 1
i += 1
OUT = Border.Edges
3 Likes
Mark, thank you for the info and translation attempt! I’ll take a closer look at that when I get a chance and update here with any progress.
I cleaned up the code a bit:
import clr
clr.AddReference('ProtoGeometry')
from Autodesk.DesignScript.Geometry import Line
from math import sqrt, hypot
points, alpha = IN
OUT = []
pLen = len(points)
alpha2 = alpha * alpha
if pLen < 2:
raise Exception('AlphaShape needs at least 2 points')
for i, p in enumerate(points):
for j in xrange(i+1, pLen):
if p.IsAlmostEqualTo(points[j]):
raise Exception('AlphaShape needs pairwise distinct points')
dist = hypot(p.X - points[j].X, p.Y - points[j].Y)
if (dist > 2 * alpha) : continue #circle fits between points ==> p_i, p_j can't be alpha-exposed
x1, y1, x2, y2 = p.X, p.Y, points[j].X, points[j].Y
midX, midY = (x1 + x2) / 2, (y1 + y2) / 2
#find two circles that contain p_i and p_j; note that center1 == center2 if dist == 2*alpha
alphaDist = sqrt(alpha2 - (dist / 2) ** 2)
deltaX, deltaY = (x2 - x1) / dist, (y1 - y2) / dist
c1x, c1y = midX + alphaDist * deltaY, midY + alphaDist * deltaX
c2x, c2y = midX - alphaDist * deltaY, midY - alphaDist * deltaX
#check if one of the circles is alpha-exposed, i.e. no other point lies in it
c1_empty = True
c2_empty = True
for k in xrange(pLen):
if i == k or j == k:
continue
if ((c1x - points[k].X) * (c1x - points[k].X) + (c1y - points[k].Y) * (c1y - points[k].Y) < alpha2):
c1_empty = False
if ((c2x - points[k].X) * (c2x - points[k].X) + (c2y - points[k].Y) * (c2y - points[k].Y) < alpha2):
c2_empty = False
if not c1_empty and not c2_empty:
break
if c1_empty or c2_empty:
OUT.append(Line.ByStartPointEndPoint(p, points[j]))
Another interesting question is, how do you determine the alpha (scoop) size? It’s correlated to the distance to the nearest neighbors but can’t be deducted exactly from it:
26 Likes
Fantastic, thanks 
1 Like
Dimitar,
First of all thank you for cleaning Mark’s script and sharing yours! This is exactly what I was looking for.
Yes, determining the alpha size is still an issue. Especially in applications with a large amount of points where testing different alpha values takes a long time to generate.
The reason I’m working on this script is an attempt to save time from having to manually draw lines in CAD; however, if the script for creating a concave hull requires fine tuning of the alpha value it could take more time to do so with the script than just manually doing it.
Do you think it is possible to find a solution that automatically finds the value (or maybe a range to pick from) of the alpha to get the desired result? I see you attempted that, but unfortunately it did not work in my application.
Also, do you think it is possible that another method of this algorithm may be more accurate? In my research I saw there was different methods (e.g. closest neighbor approach) to achieve a concave hull.
Would love to hear your thoughts and perhaps find a working solution! Thanks again.
1 Like
Ha, @Mark.Ackerley, uni + work does leave for very little time, but I was intrigued.
I know this has been answered, but here is another method that uses the Graham Scan method for “Gift Wrapping”. The best method is Chans algorithm from what I’ve read though as it is a combination of the Jarvis March and Graham Scan methods.
Here is the Python Code ported over from Thomas Switzers GitHub (I take zero credit for this)…
import clr
clr.AddReference('ProtoGeometry')
from Autodesk.DesignScript.Geometry import *
def ToVals(pt):
return [pt.X,pt.Y,pt.Z]
def ToPts(ptVals):
return Point.ByCoordinates(ptVals[0],ptVals[1],ptVals[2])
def CreateOutline(pts):
crvs = []
i = 0
while i < len(pts):
if i == len(pts)-1:
crvs.append(Line.ByStartPointEndPoint(pts[i],pts[0]))
else:
crvs.append(Line.ByStartPointEndPoint(pts[i],pts[i+1]))
i = i+1
return PolyCurve.ByJoinedCurves(crvs)
'''
Taken from Thomas Switzers awesome GitHub repo and modified to suit Dynamo Context
https://gist.github.com/arthur-e/5cf52962341310f438e96c1f3c3398b8
I take zero credit for this work other than porting it across.
'''
def convex_hull_graham(points):
'''
Returns points on convex hull in CCW order according to Graham's scan algorithm.
By Tom Switzer <thomas.switzer@gmail.com>.
'''
TURN_LEFT, TURN_RIGHT, TURN_NONE = (1, -1, 0)
def cmp(a, b):
return (a > b) - (a < b)
def turn(p, q, r):
return cmp((q[0] - p[0])*(r[1] - p[1]) - (r[0] - p[0])*(q[1] - p[1]), 0)
def _keep_left(hull, r):
while len(hull) > 1 and turn(hull[-2], hull[-1], r) != TURN_LEFT:
hull.pop()
if not len(hull) or hull[-1] != r:
hull.append(r)
return hull
points = [ToVals(p) for p in points]
points = sorted(points)
l = reduce(_keep_left, points, [])
u = reduce(_keep_left, reversed(points), [])
points = l.extend(u[i] for i in range(1, len(u) - 1)) or l
points = [ToPts(p) for p in points]
return CreateOutline(points)
OUT = convex_hull_graham(IN[0])
Cheers,
Dan
18 Likes
Thanks Dan great work 
Nah, scratch that… I thought you were after convex Hull… This does convex not concave… My bad 
Still cool to see, thanks for posting
1 Like
Sorry, I don’t have anything better. The only alternative is to use a while loop / goal seek type of approach, where you try to minimize the alpha size, while still producing a single continuous external loop.
1 Like
No problem. Thank you for all the help so far, I really appreciate it!
Looking forward to seeing you around the forums.
There is an issue with the above Python code (the one by Dimitar and Mark) if you use it for vertical point sets
the hypothemus (dist) will return zero if two points are vertically above eachother.
And then a division error comes.
Great topic though. Interesting stuff
Comes in handy when working with pointclouds
2 Likes
Thanks its magic…!!
1 Like
Is it possible to adjust it for tolerance in the code?