# Pathfinding on 2D regions

How can I find a shortest or near optimal route between two points where the route is constrained within a 2D region?

First, consider the following bundle of lines:

SeedRandom[1];
points = RandomPoint[Disk[], 70];
nf = Nearest[points];
lines = Line /@ Partition[points, 2];
start = First[nf[{0, -1}]];
end = First[nf[{0, 1}]];
Graphics[{lines, Blue, PointSize[Large], Point[start], Red, Point[end]}]


To solve this one could build a graph where the intersections are the vertices. However, what if we have a more complicated combination of regions like the following:

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

Disk[#1, #2] &, {RandomPoint[Disk[], numdisks],
RandomReal[1/5, numdisks]}];

Translate[#1, #2] &, {RandomPolygon[8, numpolys,
DataRange -> {-.15, .15}], RandomPoint[Disk[], numpolys]}];

Graphics[{
disks, polygons, PointSize[Large], Cyan, Point[{-.4, .9}], Magenta,
Point[{-.8, -.6}]
}]


There should be some path composed of line segments that gets us from the cyan dot to the magenta dot. I'd like to solve this particular example in an agnostic sense without considering any special properties of the underlying primitives. In other words, we're just given a single region like ImageMesh[ColorNegate[Graphics[{polygons, disks}]]] and there's no way to break it down further.

• Ideas at the moment: use Thinning[ColorNegate[ImageMesh[...]]] or SkeletonTransform or MorphologicalGraph to build a connectivity graph, then find a path from the vertex nearest the start point to the vertex nearest the end point. It's not very near optimal, and annoyingly requires converting to an image, which means I have to convert coordinates. To visualize: img = Graphics[{polygons, disks}]; ImageAdd[ColorReplace[Thinning[ColorNegate[img]], White -> Cyan], img] Sep 30 '20 at 13:49
• Perhaps relevant: community.wolfram.com/groups/-/m/t/1144481 Oct 1 '20 at 12:04
• Perhaps relevant. Oct 1 '20 at 13:30
• What you can do is to transform this into an image-based problem. With an image, you can do a flood fill which tracks the distance starting from one of the points. When you reach the second point, you can backtrack to find the exact path taken. Look e.g. on this gif of the A* algorithm. Oct 1 '20 at 15:21
• Related to the flood fill approach, I will remark that this can be treated as a maze problem. (I used a flood fill type of method in a maze path-finding blog several years ago. I did not realize the method had a name. Good grief, that blog is going to turn ten years old soon.) Oct 4 '20 at 15:43

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 functionality is now available via the resource function ResourceFunction["RegionFindShortestPath"]. We can use it in this case as (the original code of this answer can be found at the bottom):

(* generate the region *)
SeedRandom[1];
numdisks = 60;
numpolys = 40;

Disk[#1, #2] &, {RandomPoint[Disk[], numdisks],
RandomReal[1/5, numdisks]}];
translatePoly[poly_, pos_] :=
Polygon[# + pos & /@ poly[[1]], poly[[2]]];
polygons =
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]}]

(* create the mesh *)
mesh = DiscretizeRegion[RegionUnion[Join[polygons, disks]]];

(* create a RegionShortestPathFunction *)
spf = ResourceFunction["RegionFindShortestPath"][mesh]

(* use the function *)
Manipulate[
Show[
mesh,
Graphics[{Thick, Red, Dynamic@Line@spf[p1, p2]}]
],
{{p1,start}, Locator}, {{p2,end}, Locator}
]


# Explanation

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[
{{GrayLevel@0.5,
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], GrayLevel@0.7],
ElisionsDumpcommonGraphicsOptions
]

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

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

Disk[#1, #2] &, {RandomPoint[Disk[], numdisks],
RandomReal[1/5, numdisks]}];
translatePoly[poly_, pos_] :=
Polygon[# + pos & /@ poly[[1]], poly[[2]]];
polygons =
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.

• Wonderful answer and a nice intuitive way of thinking about it! After the pre-processing, the responsiveness with the locators is really good. Oct 2 '20 at 22:01
• Thanks for the accept! - in case you didn't notice yet (the timing was a bit unfortunate, since the edit was almost at the same time as your comment), I have improved the code to run twice as fast (and it should be easier to understand as well) Oct 2 '20 at 22:03
• @flinty Ok, so I looked into the idea with convex polygons a bit, and I don't think it will help much: For the region in the question, only around 1% of the lines are inside one of the convex polygons returned by PolygonDecomposition, so the speedup will be marginal at best Oct 5 '20 at 10:45
• @flinty After some tinkering, the code is now an additional 8 times faster, the whole thing (apart from region creation and compilation of the functions) takes less than 3 seconds now: It occurred to me that we only need to consider inwards pointing corners of the region, since the others will never be part of the shortest path. This drastically reduces the number of possible lines, leading to the speedup. Oct 6 '20 at 12:21
• @flinty The function is now (finally) available via the WFR, as ResourceFunction["RegionFindShortestPath"], see also the update in the answer Jun 7 at 22:11

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;

Disk[#1, #2] &, {RandomPoint[Disk[], numdisks],
RandomReal[1/5, numdisks]}];
translatePoly[poly_, pos_] :=
Polygon[# + pos & /@ poly[[1]], poly[[2]]];
polygons =
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]}]


• This gives a Splice::string message and the path is not displayed! Mathematica shows a page Splice is superseded by FileTemplate Oct 1 '20 at 11:59
• @SteffenJaeschke are you using 12.1.1.0 ? I do not get any messages from a fresh kernel. Oct 1 '20 at 12:47
• ^ yes this must be true. You are not using v12.1 as Splice took on a completely different meaning in 12.1 - see here: mathematica.stackexchange.com/questions/221414/… You can replace this line {start, Splice[getcoords /@ path], end} with Join[{start},getcoords /@ path, {end}] Oct 1 '20 at 12:54
• getcoords does not operate because AnnotationValues gives always the output Failed!. Look at my solution in my answer how to do it, please. Oct 1 '20 at 13:53

here is a first try. It can be improved, but one has to start somewhere. The following program takes a region, a start and end point and a step size. It makes a plot of the path (red) and the tried points (green). But take car, because there are no preconditions to exploit, it takes a lot of steps. It is a recursive program, therefore we need to enlarge "$RecursionLimit". Take care with "stepsize", make it as big as reasonable, otherwise the number of steps will explode. Further, the difference in x and y coordinates between start and endpoint must be a multiple of the stepsize. Otherwise the end is never found. And for simplicity, the start point should be chosen above the end point. Note also that the path can be pretty big, so it make no sense to print it. The idea is, that you will do something with it besides printing. Further, you will need some real region, not translated polygons, that will work with "RegionMember". Here is the program, have fun: getPath[region_, start_, end_, stepsize_] := Module[{path = {pos = start}, step = stepsize, wrong = {}, remem, search}, If[Mod[(end - start)/step, 1] != {0, 0}, Print["Difference between end and start must be a multiple of \ stepsize."]; Return[]]; remem = RegionMember[region]; search[pos_] := Module[{}, If[ Norm[pos - end] < 0.001, Return[{}]]; Which[ tp = pos + {0, -1} step; ( ! MemberQ[path, tp]) && ( ! MemberQ[wrong, tp]) && remem[tp], AppendTo[path, tp]; search[tp], tp = pos + {1, 0} step ; ( ! MemberQ[path, tp]) && ( ! MemberQ[wrong, tp]) && remem[tp], AppendTo[path, tp]; search[tp], tp = pos + {-1, 0} step; ( ! MemberQ[path, tp]) && ( ! MemberQ[wrong, tp]) && remem[tp], AppendTo[path, tp]; search[tp], True, AppendTo[wrong, path[[-1]]]; path = Delete[path, -1]; If[path == {}, Return[{}]]; search[path[[-1]]]; ]]; search[start]; Show[Region[region], Graphics[{Green, Point[wrong], Thick, Red, Line[path], Black, PointSize[0.03], Point[end] , Point[start]}], PlotRange -> {{-1, 4}, {-1, 5}}, Axes -> True] // Print; path ];  We create some arbitrary region, choose a start and end point and let the program search a path. region = RegionUnion[Disk[{0, 3.2}], Disk[{0.9, 2.2}, 0.5], Disk[{1.9, 3.}, 0.8], Disk[{2.5, 1.8}, 0.6], Disk[{1.8, .6}, .9], Disk[{0, 0}]]; start = {0., 3.2}; end = {0, 0}; stepsize = 1/10;$RecursionLimit = 10^4;
path = getPath[region, start, end, stepsize];


RegionUnion[Disk[{0, 0}, 2], Disk[{3, 0}, 2]];

Region[%]


For a simple circle and a point:

RegionDistance[Disk[{0, 0}, 2], {3, 0}]


1

Graphics[{Disk[{0, 0}, 2], Point[{3, 0}], Red,
Line[{{0, 0}, {3, 0}}]}]


If the main intention remains to work with Random-function the ideas of @flinty are not bad to look whether these are connected and a path exists.

This is the generated approach to a plane geometric arrangement of Circle and Polygon. For each the center is known and a Sort or else is easily done.

This process have to be repeated with care for each smaller set that works.

Dealing with Transform fails for BooleanRegion.

ℜpolygon =
Region@RegionUnion[
Table[Polygon[
Plus[cent[[i]], #] & /@ RandomReal[{-0.15, 0.15}, {8, 2}]], {i,
30}]]
ℜcircle =
Region@RegionUnion[
Disk[#1, #2] &, {RandomPoint[Disk[], numdisks],
RandomReal[1/5, numdisks]}]]

ℜcomp =
Region@RegionUnion[ℜpolygon, ℜcircle]


But RegionNearest and RegionDistance do not work, are not defined for BooleanRegion.

RegionDistance[
Region@RegionUnion[{Disk[{0, 0}, 2], Disk[{1, 1}, 2],
Disk[{1, -1}, 2]}], {3, 0}]

Graphics[{Disk[{0, 0}, 2], Disk[{1, 1}, 2], Disk[{1, -1}, 2],
Point[{3, 0}], Red,
Line[{{Sqrt[(3/2 + 1/10 (-5 - 4 Sqrt[5]))^2 + (-3 +
1/5 (5 + 4 Sqrt[5]))^2], 0}, {3, 0}}]}]


This too must be a BooleanRegion.

Mathematica V12 has the built-in RandomInstance and GeometricScene

RandomInstance[GeometricScene[{a, b, c, d, g, e, f}, {
a == {-1/2, 0}, b == {1/2, 0}, Line[{f, a, b, e}],
p0 == Polygon[{e, g, f}],
p1 == Style[Polygon[{a, c, b}], Yellow],
p2 == Style[Polygon[{b, d, c}], Magenta],
p3 == Style[Polygon[{d, c, g}], Green],
p4 == Style[Polygon[{g, c, a}], Blue],
p5 == Style[Polygon[{e, b, d}], Purple],
p6 == Style[Polygon[{g, a, f}], Orange],
GeometricAssertion[{p0, p1, p2, p3, p4, p5, p6}, "Similar"]}],
RandomSeeding -> 4]


It has the built-in GeometricAssertion with offers a generative process to construct the path along with the objects. And this allows nicer and more realistic random polygons. And it provides a description of the paths in the plane and it handles more geometric relations i. e. SimplePolygonQ.

Splice might be reintroduced on V12.1 and later. I can be found in other answers to question on this community. Sequence@@ old style.

AnnotationValue is not in my documentation of Mathematica V12.0.0. But AnnotationValue works on V12.0.0. So this is a built-in without value in this question: Failed.

So the rest does not work.

This works on V12.0.0 and alike:

Show[Graphics[{mesh, PointSize[Large], Cyan, Point[start], Magenta,
Point[end]}], Subgraph[cellGr, PathGraph[path]]]


Using

connectedCells[cells1_, cells2_] :=
Length[Intersection[cells1, cells2]] == 1


in the above use code gives:

Show[Graphics[{mesh, PointSize[Large], Cyan, Point[start], Magenta,
Point[end]}],
Subgraph[cellGr, PathGraph[path], EdgeStyle -> {Thick, Green}]]


This is faster, but still has the problems in the lower left part of the DiscretizeRegion.

mesh = DiscretizeRegion[RegionUnion[Join[polygons, disks]], MaxCellMeasure -> 1]

mesh = DiscretizeRegion[RegionUnion[Join[polygons, disks]],
MaxCellMeasure -> {"Length" -> 1/15}, PrecisionGoal -> None];
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]] == 1
cellGr = RelationGraph[connectedCells[#1, #2] &, cells,
VertexCoordinates -> meshcentroids];
path = FindShortestPath[cellGr, startcell, endcell];

Show[Graphics[{mesh, PointSize[Large], Cyan, Point[start], Magenta,
Point[end]}], Subgraph[cellGr, PathGraph[path]]]


This shows that adaptive meshing has to be replaced by a fine mesh of regular density for the most optimal path if triangularization is used. Every corner, every smaller trespassing, every extension pointing inward or outward attracts the triangulation mesh, and there induces an oscillation in the shortest paths.

A regular does lead to oscillations as long as the cell measure is high. For smaller once the shortest stops from oscillating and gets direct and stays inside the meshed blob. If the boundary would be taken better into account the shortest might stay more in the blobs center region.

It is a compromise between time and directedness how short the shortest path will be.

• Nothing in this answer helps me find a short path within the region. If you know how to do this, please show your method and the calculated path drawn on top of the region. Sep 30 '20 at 17:23
• Apart from allowing the connectedCells[cells1_, cells2_] := Length[Intersection[cells1, cells2]] == 1 which considers triangles connected if they meet at an edge, and making some adjustments so it works in 12.0, your edits are more or less a direct copy and paste of my answer with little to no improvement in path quality or performance. Oct 1 '20 at 14:14
• OK. This make use of some built-ins compose by Your post, but I make it work under V12.0.0 and more general for other versions. I selected some different built-ins that are more efficient but sorry for that not faster. But this solution is fast Oct 1 '20 at 14:31
• That's fair enough. If you remove the first half of your answer - everything up to and including GeometricAssertion` - which I don't believe is relevant, then I will retract my downvote. Oct 1 '20 at 14:34