Mathematica Asked on April 21, 2021
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;
disks = MapThread[
Disk[#1, #2] &, {RandomPoint[Disk[], numdisks],
RandomReal[1/5, numdisks]}];
polygons = MapThread[
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.
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:
insideQ
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.
Correct answer by Lukas Lang on April 21, 2021
Start with this:
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[
MapThread[
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.
Answered by Steffen Jaeschke on April 21, 2021
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];
Answered by Daniel Huber on April 21, 2021
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]}]
Answered by flinty on April 21, 2021
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP