Qual è il modo più efficace per trovare l'intersezione di una linea e un cerchio in python?

Question

Ho un poligono composto da molti punti. Voglio trovare l'intersezione del poligono e un cerchio. Fornendo il centro del cerchio di [x0, y0] e il raggio di r0, ho scritto una funzione approssimativa per risolvere semplicemente l'equazione quadratica del cerchio e una linea. Ma che dire dell'efficienza di trovare l'intersezione di ogni segmento di linea del poligono uno ad uno? Esiste un modo più efficiente?

So che sympy fornisce già la funzione per ottenere le intersezioni di diverse geometrie. Ma anche e l'efficienza della libreria esterna come sympy rispetto a calcolarla con la mia funzione, se voglio occuparmi di molti poligoni?

def LineIntersectCircle(p,lsp,lep): 
# p is the circle parameter, lsp and lep is the two end of the line 
    x0,y0,r0 = p 
    x1,y1 = lsp 
    x2,y2 = esp 
    if x1 == x2: 
    if abs(r0) >= abs(x1 - x0): 
     p1 = x1, y0 - sqrt(r0**2 - (x1-x0)**2) 
     p2 = x1, y0 + sqrt(r0**2 - (x1-x0)**2) 
     inp = [p1,p2] 
     # select the points lie on the line segment 
     inp = [p for p in inp if p[1]>=min(y1,y2) and p[1]<=max(y1,y2)] 
    else: 
     inp = [] 
    else: 
    k = (y1 - y2)/(x1 - x2) 
    b0 = y1 - k*x1 
    a = k**2 + 1 
    b = 2*k*(b0 - y0) - 2*x0 
    c = (b0 - y0)**2 + x0**2 - r0**2 
    delta = b**2 - 4*a*c 
    if delta >= 0: 
     p1x = (-b - sqrt(delta))/(2*a) 
     p2x = (-b + sqrt(delta))/(2*a) 
     p1y = k*x1 + b0 
     p2y = k*x2 + b0 
     inp = [[p1x,p1y],[p2x,p2y]] 
     # select the points lie on the line segment 
     inp = [p for p in inp if p[0]>=min(x1,x2) and p[0]<=max(x1,x2)] 
    else: 
     inp = [] 
    return inp 

Answer 1

Immagino che la tua domanda riguardi come in teoria farlo nel modo più veloce. Ma se vuoi farlo velocemente, dovresti davvero usare qualcosa che è scritto in C/C++.

Sono abbastanza abituato a Shapely, quindi fornirò un esempio di come eseguire questa operazione con questa libreria. Ci sono molte librerie di geometrie per Python. Li elencherò alla fine di questa risposta.

from shapely.geometry import LineString 
from shapely.geometry import Point 

p = Point(5,5) 
c = p.buffer(3).boundary 
l = LineString([(0,0), (10, 10)]) 
i = c.intersection(l) 

print i.geoms[0].coords[0] 
(2.8786796564403576, 2.8786796564403576) 

print i.geoms[1].coords[0] 
(7.121320343559642, 7.121320343559642) 

cosa è un po contatore po intuitivo ben fatto è che i cerchi sono i boundries di punti con aree tampone. Questo è il motivo per cui lo faccio p.buffer(3).boundry

anche l'intersezione i è una lista di forme geometriche, entrambi punti in questo caso, è per questo che devo ottenere entrambi da i.geoms[]

C'è another Stackoverflow question che va nei dettagli su queste librerie per gli interessati.

• SymPy
• CGAL Python bindings
• PyEuclid
• PythonOCC
• Geometry-simple

EDIT perchè commenti: o

Shapely si basa n GEOS (trac.osgeo.org/geos) che è costruito in C++ e considerevolmente più veloce di qualsiasi cosa tu scriva in modo nativo in python. SymPy sembra essere basato su mpmath (mpmath.org) che sembra essere anche python, ma sembra avere un sacco di matematica abbastanza complessa in esso integrata. Implementare ciò stesso potrebbe richiedere molto lavoro e probabilmente non sarà così veloce come le implementazioni di GEOS C++.

Answer 2

Penso che la formula che si utilizza per trovare le coordinate delle due intersezioni non possa essere ulteriormente ottimizzata. L'unico miglioramento (che è numericamente importante) è distinguere i due casi: |x_2-x_1| >= |y_2-y_1| e |x_2-x1| < |y_2-y1| in modo che la quantità k sia sempre compresa tra -1 e 1 (nel calcolo si possono ottenere errori numerici molto elevati se | x_2-x_1 | è molto piccolo). È possibile scambiare x-se y-s per ridurre un caso all'altro.

È anche possibile eseguire un controllo preliminare: se entrambi i punti terminali sono interni al cerchio, non vi sono intersezioni. Calcolando la distanza quadrata tra i punti e il centro del cerchio, questa diventa una formula semplice che non utilizza la funzione radice quadrata. L'altro controllo: "se la linea è fuori dal cerchio" è già implementato nel tuo codice e corrisponde al delta < 0. Se hai molti segmenti piccoli questi due controlli dovrebbero dare una risposta di scelta rapida (nessuna intersezione) nella maggior parte dei casi.

Answer 3

Un basso costo partizione spaziale potrebbe essere quella di suddividere il piano in 9 parti

Ecco un diagramma scadente. Immagina, se vuoi, che le linee stanno semplicemente toccando il cerchio.

 
    | | 
__|_|__ 
__|O|__ 
    | | 
    | | 

8 delle aree che ci interessa sono circondano il cerchio. Il quadrato al centro non è molto utile per un test economico, ma puoi inserire un quadrato di r/sqrt(2) all'interno del cerchio, quindi gli angoli toccano semplicemente il cerchio.

Consente etichettare le aree

 
A |B| C 
__|_|__ 
D_|O|_E 
    | | 
F |G| H 

e la piazza del r/sqrt(2) nel centro chiameremo J

Chiameremo l'insieme dei punti in piazza del centro mostrato nel diagramma che aren' t in J, Z

Etichetta ogni vertice del poligono con il suo codice di lettera.

Ora possiamo vedere rapidamente

 
AA => Outside 
AB => Outside 
AC => Outside 
... 
AJ => Intersects 
BJ => Intersects 
... 
JJ => Inside 

Questo può trasformato in una tabella di ricerca

Quindi, a seconda del set di dati, potrebbe essere stato salvato da soli un carico di lavoro. Qualunque cosa con un endpoint in Z dovrà comunque essere testata.