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

Hi Yi

I very much appreciated finding your 2D alpha shape demo in Mathematica. I had been trying to turn the mathematics in the paper “On the Shape of a set of Points in the Plane” by Herbert Edelsbrunner, D avid G Kirkpatrick and Raimund Seidel, IEEE IT-29, #4, 1983 into Mathematica with a lot of work and limited success. It is not clear to me that the test for whether an edge belongs to the alpha hull (for a given value of alpha) in this paper is the same as your elegantly simple test, i.e., one circle contains points and the other does not. This makes the code equally elegant and simple, whereas the Edelbrunner paper lead me to try to find “furthest point” Voronoi regions as well as the usual “nearest point” regions given by Mathematica’s Computational Geometry package. Can you tell me the source of your criterion, as I would like to cite it.

Thanks very much,

Barrie

Hi Barrie,

That demo was not intended to be a precise implementation of Edelsbrunner et al 1983, but was written to gain some insight of the alpha-shape PHILOSOPHY during my research on finite element, and I have not inquired its robustness and stability. Computation geometry is not my field :). For complete implementation, you may check out, among others:

H. Edelsbrunner and E.P. Mucke, Three-dimensional alpha shapes. 1994;

Sergio R. Idelsohn et al, Polyhedrization of an arbitrary 3D point set 2003;

GID user’s guide from CIMNE.

Good luck

Yi