Pathfinding on 2D regions

Here's an approach that should produce the globally optimal solution (code below):

After some preprocessing, the performance is real-time capable as shown in the gif. The preprocessing needs to be run once for each region, but takes less than 3 seconds on my machine for the region in the question.

The idea is that every shortest path will essentially consist of straight lines between points on the boundary of the region (and of course the start and end point). To see this, imagine being in a room with the shape of the region, and your candidate shortest path is marked out with a string: If you now pull on the string (to minimize the path length taken by the string), the string will be caught by some corners of the room, but will go in straight lines in between. At this point we also note that only corners pointing inward need to be considered: No shortest path will ever go to an outwards facing corner of the region, as can again be seen from the analogy with the string.

The implementation selects all inwards pointing corners in pointData (which also contains data for the function insideQ described below) and generates a list of all possible lines between any such points, and then selects those that are inside the region (this is the step that will take a while, since there are ~25000 lines to check for the region above). To get the actual path from start to end, we need to add all lines from those two points to any inwards pointing boundary point, but that list is way shorter and thus it can be computed in real time.

The tricky thing is to get a function that can quickly check whether a line is inside the region or not - the built-in region functionality is way too slow (and buggy) unfortunately, so we need a custom solution.

This is done by the functions lineWithinQ, intersectingQ and insideQ:

• insideQ checks whether the line under test points inwards from the edge of the boundary by essenitally computing the triple product of the two adjancent edge vectors and the line in question. We also compile the function for maximum performance.

• intersectingQ checks whether the line under test intersects with any of the boundary lines (touching the line does not count). The function effectively solves for the intersection of the two lines (given their endpoints) and verifies that the intersection is indeed between the endpoints. For maximum performance, this function is compiled and aborts as soon as an intersection is found

• Finally, lineWithinQ then checks whether a line is inside the region in two steps:

  • First, check whether the line points into the region at all with insideQ
  • Second, check whether the line crosses the boundary at any point with intersectingQ (remember that touching doesn't count)

Since the functions work only for lines between points on the border, adding the start and end point is done a bit differently (as seen by the handling of start and end inside the code of RegionShortestPathFunction below): We first filter lines from any boundary point to the start/end using lineWithinQ, since the function still works as long as the first point is on the boundary (insideQ checks whether the line points into the region only looking from the starting point of the line). To check whether the line straight from start to end is valid, we simply check whether it intersects the boundary at all.

Module[
 {cond, l, i},
 cond = Unevaluated@FullSimplify[0 < t < 1 && 0 < u < 1] /. 
   First@Solve[{t, 1 - t}.{{x1, y1}, {x2, y2}} == {u, 
        1 - u}.{{x3, y3}, {x4, y4}}, {t, u}];
 cond = cond /. 
   Thread[{x1, y1, x2, y2} -> Table[Indexed[l, {i, j}], {j, 4}]];
 cond = cond /. Thread[{x3, y3} -> Table[Indexed[p1, i], {i, 2}]];
 cond = cond /. Thread[{x4, y4} -> Table[Indexed[p2, i], {i, 2}]];
 With[
  {cond = cond},
  intersectingQ = Compile @@ Hold[
     {{l, _Real, 2}, {p1, _Real, 1}, {p2, _Real, 1}},
     Module[{ret = False}, 
      Do[If[cond, ret = True; Break[]], {i, Length@l}]; ret],
     CompilationTarget -> "C", RuntimeAttributes -> {Listable}, 
     Parallelization -> True
     ]
  ]
 ]

Module[
 {cond, x1, y1, z1, x2, y2, v1, v2},
 cond = {x1, y1, z1}.Append[Normalize@{x2, y2}, 1] > 0 /. 
    Abs -> RealAbs // FullSimplify[#, x2^2 + y2^2 > 0] &;
 cond = cond /. Thread[{x1, y1, z1} -> Table[Indexed[v1, i], {i, 3}]];
 cond = cond /. Thread[{x2, y2} -> Table[Indexed[v2, i], {i, 2}]];
 insideQ = Compile @@ {
    {{v1, _Real, 1}, {v2, _Real, 1}},
    cond,
    CompilationTarget -> "C", RuntimeAttributes -> {Listable}, 
    Parallelization -> True
    }
 ]

lineWithinQ[lineData_, {{p1_, v1_}, {p2_, _}}] :=
 insideQ[v1, p2 - p1] && ! intersectingQ[lineData, p1, p2]

Options[RegionFindShortestPath] = {"MonitorProgress" -> True};

RegionFindShortestPath[region_?MeshRegionQ, start : {_, _}, end : {_, _}, opts : OptionsPattern[]] :=
 RegionFindShortestPath[region, start, opts][end]
RegionFindShortestPath[region_?MeshRegionQ, start : {_, _}, opts : OptionsPattern[]] :=
 RegionFindShortestPath[region, opts][start]

RegionFindShortestPath[region_?MeshRegionQ, OptionsPattern[]] :=
 Module[
  {lines, lineData, pointData, pathData},
  lines = MeshPrimitives[RegionBoundary@region, 1][[All, 1]];
  lineData = Catenate /@ lines;
  pointData = Cases[(* select inwards pointing corners *)
     {p_, {__, z_} /; z > 0, c_} :> {p, c}
     ]@Catenate[
     Transpose@{
         #[[All, 2]],
         Sequence @@ Table[
           Cross[#, {-1, -1, 1} #2] & @@@
            Partition[
             Append[z]@*Normalize /@ Subtract @@@ #,
             2, 1, {1, 1}
             ],
           {z, 0, 1}
           ]
         } & /@
      FindCycle[Graph[UndirectedEdge @@@ lines], \[Infinity], All]
     ];
  pathData = With[
    {expr := 
      Select[lineWithinQ[lineData, #] &]@Subsets[pointData, {2}]},
    If[OptionValue["MonitorProgress"],
      ResourceFunction["MonitorProgress"][expr, 
       "CurrentDisplayFunction" -> None],
      expr
      ][[All, All, 1]]
    ];
  RegionShortestPathFunction[pointData, lineData, 
   Join[pathData, lines]]
  ]

RegionShortestPathFunction[data__][start : {_, _}, end : {_, _}] :=
 RegionShortestPathFunction[data][start][end]
RegionShortestPathFunction[pointData_, lineData_, pathData_][start : {_, _}] :=
 RegionShortestPathFunction[pointData, lineData, Join[
   pathData,
   Select[lineWithinQ[lineData, #] &][{#, {start, {}}} & /@ 
      pointData][[All, All, 1]]
   ], start]

RegionShortestPathFunction[pointData_, lineData_, pathData_, start_][end : {_, _}] :=
 With[
  {allLines = Join[
     pathData,
     Select[lineWithinQ[lineData, #] &][{#, {end, {}}} & /@ 
        pointData][[All, All, 1]],
     If[! intersectingQ[lineData, start, end], {{start, end}}, {}]
     ]},
  Quiet@
   Check[
    FindShortestPath[
     Graph[UndirectedEdge @@@ allLines, 
      EdgeWeight -> EuclideanDistance @@@ allLines], start, end],
    {}
    ]
  ]

summaryBoxIcon = Graphics[
  {{[email protected], 
    Polygon@{{0, 0}, {0, 1}, {1, 1}, {1, -1}, {-2, -1}, {-2, 
       1.5}, {-1, 1.5}, {-1, 0}}}, {Red, 
    Line@{{0.5, 0.5}, {0, 0}, {-1, 0}, {-1.5, 1}}}, 
   AbsolutePointSize@4, Point[{0.5, 0.5}], {Point[{-1.5, 1}]}}, 
  Background -> GrayLevel[0.93], PlotRangePadding -> Scaled[0.1], 
  FrameStyle -> Directive[Thickness[Tiny], [email protected]], 
  ElisionsDump`commonGraphicsOptions
  ]

MakeBoxes[
  f : RegionShortestPathFunction[pointData_, lineData_, pathData_, 
    start_ | PatternSequence[]], fmt_] ^:=
 BoxForm`ArrangeSummaryBox[
  RegionShortestPathFunction,
  f,
  summaryBoxIcon,
  {
   BoxForm`SummaryItem@{"Corner points: ", Length@lineData},
   BoxForm`SummaryItem@{"Start set: ", Length@{start} > 0}
   },
  {
   BoxForm`SummaryItem@{"Possible segments: ", Length@pathData}
   },
  fmt
  ]

SeedRandom[1];
numdisks = 60;
numpolys = 40;

disks = MapThread[
   Disk[#1, #2] &, {RandomPoint[Disk[], numdisks], 
    RandomReal[1/5, numdisks]}];
translatePoly[poly_, pos_] := 
  Polygon[# + pos & /@ poly[[1]], poly[[2]]];
polygons = 
  MapThread[
   translatePoly[#1, #2] &, {RandomPolygon[8, numpolys, 
     DataRange -> {-.15, .15}], RandomPoint[Disk[], numpolys]}];
start = {-.4, .9};
end = {-.8, -.6};
Graphics[{disks, polygons, PointSize[Large], Cyan, Point[start], 
  Magenta, Point[end]}]
mesh = DiscretizeRegion[RegionUnion[Join[polygons, disks]]];

spf = RegionFindShortestPath[mesh]

Manipulate[
 Show[
   mesh,
   Graphics[{Thick, Red, Dynamic@Line@spf[p1, p2]}]
  ],
 {p1, Locator},
 {p2, Locator}
 ]

As demonstrated, the function can be used as RegionFindShortestPath[mesh][start,end] (where RegionFindShortestPath[mesh] gives a RegionShortestPathFunction with the precomputed information cached inside). All combinations such as RegionFindShortestPath[mesh,start,end] and RegionFindShortestPath[mesh,start][end] work as well, with as much information as possible being cached.

I came up with an unconventional and inefficient solution that may be susceptible to meshing problems and sensitive to mesh cell size, but I believe it produces a reasonably good short path. Maybe others could improve this or suggest alternative solutions:

SeedRandom[1];
numdisks = 60;
numpolys = 40;

disks = MapThread[
   Disk[#1, #2] &, {RandomPoint[Disk[], numdisks], 
    RandomReal[1/5, numdisks]}];
translatePoly[poly_, pos_] := 
  Polygon[# + pos & /@ poly[[1]], poly[[2]]];
polygons = 
  MapThread[
   translatePoly[#1, #2] &, {RandomPolygon[8, numpolys, 
     DataRange -> {-.15, .15}], RandomPoint[Disk[], numpolys]}];
start = {-.4, .9};
end = {-.8, -.6};
Graphics[{disks, polygons, PointSize[Large], Cyan, Point[start], 
  Magenta, Point[end]}]
mesh = DiscretizeRegion[RegionUnion[Join[polygons, disks]]];

cells = MeshCells[mesh, 2][[All, 1]];
prims = MeshPrimitives[mesh, 2];
meshcentroids = RegionCentroid /@ prims;
nprim = Nearest[meshcentroids -> "Index"];
startcell = cells[[First[nprim[start]]]];
endcell = cells[[First[nprim[end]]]];
connectedCells[cells1_, cells2_] := 
 Length[Intersection[cells1, cells2]] == 2
cellGr = RelationGraph[connectedCells[#1, #2] &, cells, 
   VertexCoordinates -> meshcentroids];
path = FindShortestPath[cellGr, startcell, endcell];

Show[Graphics[
  {EdgeForm[LightRed], FaceForm[LightYellow], mesh, PointSize[Large], 
   Cyan, Point[start], Magenta, Point[end]}
  ], HighlightGraph[cellGr, PathGraph[path]]
 ]

The above code finds a path in the mesh connectivity graph. That is the graph of adjacent triangles (sharing an edge) in the discretized mesh. This path is obviously very squiggly, so the following code tries to find the longest 'leaps' along the path that can skip out vertices but stay within the region:

(** from the currentPoint, try to draw a line that to the furthest 
  possible point on the path that stays within the region **)
getcoords[cell_] := AnnotationValue[{cellGr, cell}, VertexCoordinates]
pathcoords = Join[{start},getcoords /@ path, {end}];
maxiline[currentPoint_, coords_] := 
 SelectFirst[Reverse[coords], 
  Quiet[Check[RegionWithin[mesh, Line[{currentPoint, #}]], False]] &]
lpath = NestWhileList[maxiline[#, pathcoords] &, start, # != end &];
Graphics[{mesh, Red, Line[lpath], PointSize[Large], Cyan, 
  Point[start], Magenta, Point[end]}]