2D alpha shape demo in Mathematica

by Yi Zhang

This is the source code for a past post.

<< ComputationalGeometry`;

(*Loading test points*)
p = RandomReal[{0, 10}, {1000, 2}];
p = Select[p, Norm[# - {5, 5}] < 5 &]; p = Select[p, Norm[# - {7.5, 5}] > 2.5 &];
ListPlot[p, PlotRange -> Full, AspectRatio -> 1]
tri = DelaunayTriangulation[p];
tri = Thread[List @@ #] & /@ tri // Flatten[#, 1] &;
tri = Union[Sort /@ tri];

(*Main function def*)
(*0 for interior bar,1 for boundary bar, 2 for bar to be erased*)
ClearAll[f];
emptyQ[center_, alpha_, plist_] := Module[
   {empty = True, n = 1},
   While[empty && n <= Length[plist], empty = Norm[plist[[n]] - center] > alpha; n++];
   empty
   ];
f[alpha_, plist_, {id1_, id2_}] := Module[
   {p1 = plist[[id1]], p2 = plist[[id2]], center1, center2, lhalf,
    center1emptyQ, center2emptyQ, property},
   lhalf = Norm[p2 - p1]/2;
   center1 = (p2 + p1)/2 +
     Sqrt[(alpha/lhalf)^2 - 1] {{0, -1}, {1, 0}}.((p2 - p1)/2);
   center2 = (p2 + p1)/2 +
     Sqrt[(alpha/lhalf)^2 - 1] {{0, 1}, {-1, 0}}.((p2 - p1)/2);
   center1emptyQ =
    emptyQ[center1, alpha, Delete[plist, {{id1}, {id2}}]];
   center2emptyQ =
    emptyQ[center2, alpha, Delete[plist, {{id1}, {id2}}]];
   Which[center1emptyQ \[And] center2emptyQ,
    2, (! center1emptyQ) \[And] (! center2emptyQ),
    0, (center1emptyQ \[And] (! center2emptyQ)) \[Or] ((!
         center1emptyQ) \[And] center2emptyQ), 1]
   (*Print[property];
   Show[Graphics[{Red,Circle[center1,alpha],Circle[center2,alpha]}],
   ListPlot[plist,AspectRatio->1]]*)
   ] /; alpha > Norm[plist[[id1]] - plist[[id2]]]/2
f[alpha_, plist_, {id1_, id2_}] :=
 2 /; alpha < Norm[plist[[id1]] - plist[[id2]]]/2

(*plot a test result*)
test = Extract[tri, Position[f[0.33, p, #] & /@ tri, 1]];
Graphics[GraphicsComplex[p, Line[test]]]
Advertisements