DelaunayMesh in a specified closed region - creating a concave hull from a set of points

Question

I have a list of discrete points wich I want to use as nodes for creating a 2D mesh. I used DelaunayMesh and it works fine. The problem that I have is that some elements/polygons are outside of the closed region that points on the boundary create(outer red dots in the picture).

Is there a way to use DelaunayMesh in a defined region or delete polygons that are outside afterwards?

points I used in this example and my code:

data   = Import["C:\\path\\example.txt", {"Data",  {All}, {1, 2, 3, 4, 5}}];
data2D = data[[All,3;;4]];
R=DelaunayMesh[data2D];
HighlightMesh[%,Style[0, Red]]

Answer 1

What you want is the 2D alpha shape to try to get close to the outline you seek. Of course, it's no longer a Delaunay triangulation since you're deleting certain polygons from the DelaunayMesh. We'll adopt my answer from this post. Here it is for a 2D point set:

alphaShapes2D[points_, crit_] := 
 Module[{alphacriteria, del = Quiet @ DelaunayMesh @ points, tetras, 
   tetcoords, tetradii, selectExternalFaces}, 
  alphacriteria[tetrahedra_, radii_, rmax_] := 
           Pick[tetrahedra, UnitStep @ Subtract[rmax, radii], 1]; 
  selectExternalFaces[facets_] := MeshRegion[points, facets]; 
  If[Head[del] === EmptyRegion, del, tetras = MeshCells[del, 2]; 
   tetcoords = MeshPrimitives[del, 2][[All, 1]]; 
   tetradii = Quiet @ Thread[Circumsphere[tetcoords]][[All, 2]]; 
   selectExternalFaces @ alphacriteria[tetras, tetradii, crit]]]

Let's see how well it performs:

reg = alphaShapes2D[data2D, 1.1];
MeshRegion[reg, MeshCellStyle -> {{0, All} -> Red}]

Looks better but doesn't get rid of some stubborn triangles. Maybe there are better ways to do this but I don't know about them.

Update

I am just updating the code to boost it's performance. (The code is already fast but it never hurts to squeeze out some more speed). The idea is to replace the built-in CircumSphere with a Listable and Compiled version, we'll call it circumRadius2D. Here is the code:

circumRadius2D = Compile[{{v, _Real, 2}},
  With[{a = Norm[v[[1]] - v[[2]]], b = Norm[v[[1]] - v[[3]]], c = Norm[v[[2]] - v[[3]]]},
   (a b c) / Sqrt[(a + b + c) (b + c - a) (c + a - b) (a + b - c)]
   ], RuntimeOptions -> "Speed", RuntimeAttributes -> {Listable}, Parallelization -> True
  ]

And here is the updated concave hull (alpha shapes) code:

alphaShapes2DC[points_, crit_] := 
 Module[{alphacriteria, del = Quiet @ DelaunayMesh @ points, tris, 
   tricoords, triradii, getExternalFaces}, 
  alphacriteria[triangle_, radius_, rmax_] := 
           Pick[triangle, UnitStep @ Subtract[rmax, radius], 1]; 
  getExternalFaces[facets_] := MeshRegion[points, facets]; 
  If[Head[del] === EmptyRegion, del, tris = MeshCells[del, 2]; 
   tricoords = MeshPrimitives[del, 2][[All, 1]]; 
   triradii = circumRadius2D @ tricoords; 
   getExternalFaces @ alphacriteria[tris, triradii, crit]]]

This updated code is more than 3 times faster than the old version. Note that I have not compiled to C. Let's put it to test with Jason's data:

reg = alphaShapes2DC[pnts, .33];
MeshRegion[reg, MeshCellStyle -> {{0, All} -> Red}]

Answer 2

I was writing this code because when I googled "Concave Hull Mathematica" I didn't find this post. I went through the trouble of taking this code and improving the performance, and bringing it up to date with newer functions.

Then before posting it, I saw RunnyKine's code, which is definitely faster. I thought I'd go ahead and post this here though.

The point is that for every pair of points p1 and p2 in the convex hull mesh, we define two circles of radius alpha, which intersect p1 and p2. If the center of one, but not both, of these circles lie a distance greater than alpha from each point in the data, then this line is a boundary line. At the end, we create a BoundaryMeshRegion,

concaveHullRegion[points_, alpha_] := 
 Module[{dtri, outsideregion, boundaryLineQ},
  dtri = Union[
    Sort /@ Flatten[List @@@ MeshCells[DelaunayMesh[points], 1], 1]];
  outsideregion[center_, plist2_] := 
   Module[{empty = True, n = 1, 
     plist3 = SortBy[plist2, Norm[# - center] &]},
    Norm[plist3[[1]] - center] > alpha];
  boundaryLineQ[plist_, {id1_, id2_}] := Module[{p1 = plist[[id1]],
     p2 = plist[[id2]],
     center1, center2, lhalf},
    lhalf = Norm[p2 - p1]/2;
    If[lhalf > alpha,
     False,
     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);
     Xor @@ (outsideregion[#, 
          Delete[plist, {{id1}, {id2}}]] & /@ {center1, center2})
     ]
    ];
  BoundaryMeshRegion[points, 
   Line@Select[dtri, boundaryLineQ[points, #] &]]
  ]

I don't have the data that was used in the OP, but I'll use some other data

cup = << "http://pastebin.com/raw/DRSCc4ir";
pnts = cup[[All, ;; 2]] // DeleteDuplicates;
conc = concaveHullRegion[pnts, .33];
conv = DelaunayMesh[pnts];    
Show[#, ListPlot[pnts, PlotStyle -> Red], 
   ImageSize -> 500] & /@ {conc, conv}