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:
