Pathfinding su regioni 2D
Come posso trovare un percorso più breve o quasi ottimale tra due punti in cui il percorso è vincolato all'interno di una regione 2D?
Innanzitutto, considera il seguente fascio di linee:
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]}]
Per risolvere questo problema si potrebbe costruire un grafo in cui le intersezioni sono i vertici. Tuttavia, cosa succede se abbiamo una combinazione più complicata di regioni come la seguente:
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}]
}]
Dovrebbe esserci un percorso composto da segmenti di linea che ci porta dal punto ciano al punto magenta. Vorrei risolvere questo particolare esempio in senso agnostico senza considerare alcuna proprietà speciale delle primitive sottostanti. In altre parole, ci viene data solo una singola regione simile ImageMesh[ColorNegate[Graphics[{polygons, disks}]]]e non c'è modo di scomporla ulteriormente.
Risposte
Ecco un approccio che dovrebbe produrre la soluzione ottimale a livello globale (codice sotto):
Dopo un po 'di pre-elaborazione, la performance è in grado di funzionare in tempo reale, come mostrato nella gif. La preelaborazione deve essere eseguita una volta per ogni regione, ma sulla mia macchina sono necessari meno di 3 secondi per la regione in questione.
L'idea è che ogni percorso più breve sarà essenzialmente costituito da linee rette tra i punti sul confine della regione (e ovviamente il punto iniziale e finale). Per vederlo, immagina di essere in una stanza con la forma della regione e il tuo percorso più breve candidato è segnato con una stringa: se ora tiri la corda (per ridurre al minimo la lunghezza del percorso presa dalla stringa), la stringa essere catturato da alcuni angoli della stanza, ma andrà in linea retta nel mezzo. A questo punto notiamo anche che devono essere considerati solo gli angoli che puntano verso l'interno: nessun percorso più breve andrà mai verso un angolo rivolto verso l'esterno della regione, come si può ancora vedere dall'analogia con la corda.
L'implementazione seleziona tutti gli angoli che puntano verso l'interno pointData(che contiene anche i dati per la funzione insideQdescritta di seguito) e genera un elenco di tutte le possibili linee tra tali punti, quindi seleziona quelle che si trovano all'interno della regione (questo è il passaggio che richiederà un mentre, poiché ci sono ~ 25000 linee per verificare la regione sopra). Per ottenere il percorso effettivo dall'inizio alla fine, dobbiamo aggiungere tutte le linee da quei due punti a qualsiasi punto di confine che punta verso l'interno, ma quell'elenco è molto più breve e quindi può essere calcolato in tempo reale.
La cosa complicata è ottenere una funzione che possa controllare rapidamente se una linea è all'interno della regione o meno: la funzionalità della regione incorporata è purtroppo troppo lenta (e buggata), quindi abbiamo bisogno di una soluzione personalizzata.
Questo viene fatto dalle funzioni lineWithinQ, intersectingQe insideQ:
insideQcontrolla se la linea sottoposta a test punta verso l'interno dal bordo del confine calcolando essenzialmente il triplo prodotto dei due vettori adiacenti del bordo e la linea in questione. Compiliamo anche la funzione per le massime prestazioni.intersectingQcontrolla se la linea in prova si interseca con una qualsiasi delle linee di confine (toccare la linea non conta). La funzione risolve efficacemente l'intersezione delle due linee (dati i loro punti finali) e verifica che l'intersezione sia effettivamente tra i punti finali. Per le massime prestazioni, questa funzione viene compilata e si interrompe non appena viene trovato un incrocioInfine,
lineWithinQcontrolla se una linea è all'interno della regione in due passaggi:- Innanzitutto, controlla se la linea punta nella regione con
insideQ - Secondo, controlla se la linea attraversa il confine in qualsiasi punto con
intersectingQ(ricorda che toccare non conta)
- Innanzitutto, controlla se la linea punta nella regione con
Poiché le funzioni funzionano solo per le linee tra i punti sul confine, l'aggiunta del punto iniziale e finale viene eseguita in modo leggermente diverso (come si vede dalla gestione starte endall'interno del codice di RegionShortestPathFunctionseguito): Per prima cosa filtriamo le linee da qualsiasi punto di confine all'inizio / end using lineWithinQ, poiché la funzione funziona ancora fintanto che il primo punto è sul confine ( insideQcontrolla se la linea punta nella regione guardando solo dal punto iniziale della linea). Per verificare se la linea retta dall'inizio alla fine è valida, controlliamo semplicemente se interseca il confine.
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}
]
Come dimostrato, la funzione può essere utilizzata come RegionFindShortestPath[mesh][start,end](dove RegionFindShortestPath[mesh]fornisce una RegionShortestPathFunctioncon le informazioni precalcolate memorizzate nella cache). Tutte le combinazioni come RegionFindShortestPath[mesh,start,end]e RegionFindShortestPath[mesh,start][end]funzionano allo stesso modo, con quante più informazioni possibili vengono memorizzate nella cache.
Ho trovato una soluzione non convenzionale e inefficiente che potrebbe essere suscettibile a problemi di mesh e sensibile alla dimensione delle celle di mesh, ma credo che produca un percorso breve ragionevolmente buono. Forse altri potrebbero migliorare questo o suggerire soluzioni alternative:
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]]
]
Il codice precedente trova un percorso nel grafico della connettività mesh. Questo è il grafico dei triangoli adiacenti (che condividono un bordo) nella mesh discretizzata. Questo percorso è ovviamente molto ondulato, quindi il codice seguente cerca di trovare i `` salti '' più lunghi lungo il percorso che possono saltare i vertici ma rimanere all'interno della regione:
(** 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]}]
ecco un primo tentativo. Può essere migliorato, ma bisogna iniziare da qualche parte. Il seguente programma prende una regione, un punto iniziale e finale e una dimensione del passo. Crea un grafico del percorso (rosso) e dei punti provati (verde). Ma prendi la macchina, perché non ci sono presupposti da sfruttare, ci vogliono tanti passaggi. È un programma ricorsivo, quindi dobbiamo ingrandire "$ RecursionLimit". Fai attenzione a "stepsize", rendilo grande quanto ragionevole, altrimenti il numero di passaggi esploderà. Inoltre, la differenza nelle coordinate xey tra il punto iniziale e quello finale deve essere un multiplo della dimensione del gradino. Altrimenti la fine non si trova mai. E per semplicità, il punto iniziale dovrebbe essere scelto sopra il punto finale. Nota anche che il percorso può essere abbastanza grande, quindi non ha senso stamparlo. L'idea è che farai qualcosa con esso oltre alla stampa.Inoltre, avrai bisogno di alcune regioni reali, non poligoni tradotti, che funzionino con "RegionMember". Ecco il programma, divertiti:
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
];
Creiamo una regione arbitraria, scegliamo un punto di inizio e di fine e lasciamo che il programma cerchi un percorso.
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];
Inizia con questo:
RegionUnion[Disk[{0, 0}, 2], Disk[{3, 0}, 2]];
Region[%]
Per un semplice cerchio e un punto:
RegionDistance[Disk[{0, 0}, 2], {3, 0}]
1
Graphics[{Disk[{0, 0}, 2], Point[{3, 0}], Red,
Line[{{0, 0}, {3, 0}}]}]
Se l'intenzione principale rimane quella di lavorare con Random-function, le idee di @flinty non sono male per vedere se queste sono connesse e se esiste un percorso.
Questo è l'approccio generato a una disposizione geometrica piana di Circlee Polygon. Per ognuno il centro è noto e un Sortaltro è facilmente realizzabile.
Questo processo deve essere ripetuto con cura per ogni set più piccolo che funziona.
Trattare con Transformfallisce per 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]
Ma RegionNearest e RegionDistance non funzionano, non sono definiti per 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}}]}]
Anche questo deve essere un file BooleanRegion.
Mathematica V12 ha il built-in RandomInstanceeGeometricScene
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]
Ha il built-in GeometricAssertioncon offre un processo generativo per costruire il percorso insieme agli oggetti. E questo consente poligoni casuali più belli e realistici. E fornisce una descrizione dei percorsi nel piano e gestisce più relazioni geometriche es SimplePolygonQ.
Splice potrebbe essere reintrodotto su V12.1 e versioni successive. Posso essere trovato in altre risposte a domande su questa comunità. Sequenza @@ vecchio stile.
AnnotationValuenon è nella mia documentazione di Mathematica V12.0.0. Ma AnnotationValuefunziona su V12.0.0. Quindi questo è un built-in senza valore in questa domanda: Fallito.
Quindi il resto non funziona.
Funziona su V12.0.0 e simili:
Show[Graphics[{mesh, PointSize[Large], Cyan, Point[start], Magenta,
Point[end]}], Subgraph[cellGr, PathGraph[path]]]
Utilizzando
connectedCells[cells1_, cells2_] :=
Length[Intersection[cells1, cells2]] == 1
nel codice di utilizzo sopra riportato fornisce:
Show[Graphics[{mesh, PointSize[Large], Cyan, Point[start], Magenta,
Point[end]}],
Subgraph[cellGr, PathGraph[path], EdgeStyle -> {Thick, Green}]]
Questo è più veloce, ma ha ancora i problemi nella parte inferiore sinistra del file 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]]]
Ciò mostra che la mesh adattativa deve essere sostituita da una mesh fine di densità regolare per il percorso più ottimale se viene utilizzata la triangolarizzazione. Ogni angolo, ogni più piccolo sconfinamento, ogni estensione che punta verso l'interno o verso l'esterno attrae la maglia di triangolazione e lì induce un'oscillazione nei percorsi più brevi.
Un normale porta a oscillazioni fintanto che la misura della cella è alta. Per una volta più piccolo, il più breve smette di oscillare e diventa diretto e rimane all'interno del blob mesh. Se il confine fosse preso meglio in considerazione, il più corto potrebbe rimanere di più nella regione centrale dei blob.
È un compromesso tra tempo e direzione quanto sarà breve il percorso più breve.