Finding minimal region with padding and linear map to circle

Question

Suppose I have n ordered points. I would like to find a simple region that contains all the points and is close to the smallest region that contains the points except it is smooth and has some "padding"(as if the points were disks with some radius).

One way to to think about it is if the points were fixed "marbles" and one were to use shrink wrap around them all. It would sort of be the type of region I'm looking for(although I sort of want it to be "locally" convex).

Another way to think about it is a sort of concave hull but with some control parameter that changes how convex it is locally. Take two disks A and B some distance apart. The convexity parameter would control how much the region "dipped" in the middle between the points. On one hand it might look like an ellipse and on the other it might look like a figure 8.

The ability to control "how much" is what I'm looking for

I would then like to have a linear map from the circle(or convex hull) to this curve so I can "morph" the circle into this region.

For example, with such an algorithm I would create a slider and on one end it will show the convex hull or circle and as the slider moves to the other side the region starts to "shrink" around the points reducing the area but always containing the "disks" and staying rather smooth and simple.

Since the convex hull should contain this "minimal" region I'm looking for it should be easy to morph between them with just linear interpolation and still have a valid region(not minimal but smooth and containing the disks).

Answer 1

As @whuber pointed it out in a comment, you are looking for α-shapes, where the parameter α defines the "tightness" of the wrapping of the boundary of a shape. I've found an implementation by YZ here from 2009, which I've fine-tuned and made faster by compiling down frequently used functions. Slowest part is the Delaunay triangulation, but there are known methods to boost its performance.

About the padding and linear mapping to a circle: it is left for someone else to do.

emptyQ = Compile[{{center, _Real, 1}, {alpha, _Real}, {plist, _Real, 
     2}}, Module[{empty = True, n = 1, x, y},
    While[empty && n <= Length@plist, {x, y} = plist[[n]] - center; 
     empty = Sqrt[Abs[x]^2 + Abs[y]^2] > alpha; n++];
    If[empty, 1, 0]],
   Parallelization -> True,
   CompilationTarget -> "C"];

center = Compile[{{alpha, _Real}, {p1, _Real, 1}, {p2, _Real, 1}}, 
   Module[{x, y, lhalf},
    {x, y} = p2 - p1;
    lhalf = Sqrt[Abs[x]^2 + Abs[y]^2]/2;
    {(p2 + p1)/2 + Sqrt[(alpha/lhalf)^2 - 1] {{0, -1}, {1, 0}}.((p2 - p1)/2),
     (p2 + p1)/2 + Sqrt[(alpha/lhalf)^2 - 1] {{0, 1}, {-1, 0}}.((p2 - p1)/2)}],
   Parallelization -> True ,
   CompilationTarget -> "C"];

f[alpha_, plist_, {id1_, id2_}] := Module[
   {p1 = plist[[id1]], p2 = plist[[id2]], c1, c2, lhalf, c1emptyQ, c2emptyQ},
   If[alpha <= Norm[p1 - p2]/2, False,
    {c1, c2} = center[N@alpha, p1, p2];
    Sow[{c1, c2}];
    c1emptyQ = emptyQ[c1, N@alpha, Delete[plist, {{id1}, {id2}}]];
    c2emptyQ = emptyQ[c2, N@alpha, Delete[plist, {{id1}, {id2}}]];
    (c1emptyQ == 1 && c2emptyQ == 0) || (c1emptyQ == 0 && c2emptyQ == 1)
    ]];

Now do the triangulation and check for each pair if the function f returns True or False:

n = 1000; (* number of points *)
alpha = .8;

Needs@"ComputationalGeometry`";
points = Select[RandomReal[{0, 10}, {n, 2}],
    (Norm[# - {5, 5}] < 5 \[And] Norm[# - {7.5, 5}] > 2.5) &];

triangulation = 
  Union[Sort /@ 
    Flatten[Thread[List @@ #] & /@ DelaunayTriangulation@points, 1]];
centerList = Last@Last@Reap[
     boundary = Select[triangulation, f[alpha, points, #] &];
     ];

{
 ListPlot[points, PlotRange -> Full, AspectRatio -> 1, 
   ImageSize -> 300, Frame -> False, Axes -> False],
 Show[
  Graphics[{Red, [email protected], 
    Map[Circle[#, alpha] &, centerList, {2}]}, ImageSize -> 300],
  Graphics@GraphicsComplex[points, Line@boundary]
  ],
 Graphics[GraphicsComplex[points, Line@boundary], ImageSize -> 300]
 }

Table[
 boundary = Select[triangulation, f[alpha, points, #] &];
 Graphics[GraphicsComplex[points, Line@boundary], PlotLabel -> Row@{"α = ", alpha}],
 {alpha, {.1, .2, .3, .4, .6, 1.}}]