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!

# Concave Hull (Alpha Shape) from list of points with Python?

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

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 https://provingground.io/2017/10/03/lunchboxml-for-dynamo/

by @Nathan_Miller . Right now, thatâ€™s the furthest you can get.

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
```

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:

Fantastic, thanks

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.

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

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

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.

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