Ricerca di punti di intersezione di due ellissi (Python)

Question

Sto scrivendo una libreria di forme 2D di base in Python (principalmente per manipolare i disegni SVG) e non riesco a calcolare in modo efficiente i punti di intersezione di due ellissi.

Ogni ellisse è definita dalle seguenti variabili (tutti float):

c: center point (x, y) 
hradius: "horizontal" radius 
vradius: "vertical" radius 
phi: rotation from coordinate system's x-axis to ellipse's horizontal axis 

ignorando quando le ellissi sono identici, ci potrebbe essere da 0 a 4 punti di intersezione (intersezione, tangenti, parzialmente sovrapposti, parzialmente sovrapposti e internamente tangente e completamente sovrapposto).

ho trovato alcune soluzioni possibili:

• SymPy geometry module - Questo fondamentalmente solo tappi le equazioni ellisse in risolutore di SymPy. Non sono sicuro se questo abbia senso senza avere già il risolutore. (Per inciso, avrei usato SymPy invece di girare il mio, ma si comporta in modo orribile quando si tratta di galleggianti pazzi)
• How to detect if an ellipse intersects(collides with) a circle - Questo potrebbe essere probabilmente adattato per due ellissi, ma sono un po 'sfocato su come girarlo in un codice sensato.
• How Ellipse to Ellipse intersection? - La libreria dei riferimenti di risposta (CADEMIA) potrebbe avere un buon algoritmo, ma non riesco nemmeno a capire se è open source.
• Wikipedia: Intersecting Two Conics - Non ho abbastanza una comprensione di algebra lineare per comprendere questa soluzione.

Qualche suggerimento su come dovrei procedere per calcolare le intersezioni? Velocità (potrebbe essere necessario calcolare molte intersezioni) e l'eleganza è il criterio principale. Il codice sarebbe fantastico, ma anche una buona direzione per entrare sarebbe utile.

Answer 1

In matematica, è necessario esprimere le ellissi come equazioni quadratiche bivariate e risolverlo. Ho trovato un doucument. Tutti i calcoli sono nel documento, ma potrebbe richiedere del tempo per implementarlo in Python.

Così un altro metodo è quello di ravvicinare le ellissi come polilinee, e utilizzare formosa per trovare le intersezioni, ecco il codice:

import numpy as np 
from shapely.geometry.polygon import LinearRing 

def ellipse_polyline(ellipses, n=100): 
    t = linspace(0, 2*np.pi, n, endpoint=False) 
    st = np.sin(t) 
    ct = np.cos(t) 
    result = [] 
    for x0, y0, a, b, angle in ellipses: 
     angle = np.deg2rad(angle) 
     sa = np.sin(angle) 
     ca = np.cos(angle) 
     p = np.empty((n, 2)) 
     p[:, 0] = x0 + a * ca * ct - b * sa * st 
     p[:, 1] = y0 + a * sa * ct + b * ca * st 
     result.append(p) 
    return result 

def intersections(a, b): 
    ea = LinearRing(a) 
    eb = LinearRing(b) 
    mp = ea.intersection(eb) 

    x = [p.x for p in mp] 
    y = [p.y for p in mp] 
    return x, y 

ellipses = [(1, 1, 2, 1, 45), (2, 0.5, 5, 1.5, -30)] 
a, b = ellipse_polyline(ellipses) 
x, y = intersections(a, b) 
plot(x, y, "o") 
plot(a[:,0], a[:,1]) 
plot(b[:,0], b[:,1]) 

e l'uscita:

Ci vuole circa 1,5ms sul mio PC.

Answer 2

guardando sympy Penso che abbia tutto il necessario. (Ho cercato di offrire esempi di codici, purtroppo, non sono riuscito ad installare sympy con gmpy2 invece di pitone inutili built-in matematica)

si ha:

• un ellipse con rotate metodo, che può essere intersected con altre ellissi
• o una funzione arbitrary intersection che richiede un numero variabile di .. che chiamano "entità geometriche".

dai loro esempi, penso che sia sicuramente possibile intersecare due ellissi:

>>> from sympy import Ellipse, Point, Line, sqrt 
>>> e = Ellipse(Point(0, 0), 5, 7) 
... 
>>> e.intersection(Ellipse(Point(1, 0), 4, 3)) 
[Point(0, -3*sqrt(15)/4), Point(0, 3*sqrt(15)/4)] 
>>> e.intersection(Ellipse(Point(5, 0), 4, 3)) 
[Point(2, -3*sqrt(7)/4), Point(2, 3*sqrt(7)/4)] 
>>> e.intersection(Ellipse(Point(100500, 0), 4, 3)) 
[] 
>>> e.intersection(Ellipse(Point(0, 0), 3, 4)) 
[Point(-363/175, -48*sqrt(111)/175), Point(-363/175, 48*sqrt(111)/175), 
Point(3, 0)] 
>>> e.intersection(Ellipse(Point(-1, 0), 3, 4)) 
[Point(-17/5, -12/5), Point(-17/5, 12/5), Point(7/5, -12/5), 
Point(7/5, 12/5)] 

edit: dal ellisse generale (ax^2 + bx + cy^2 + dy + exy + f) non è supportato in sympy,

si dovrebbero costruire equazioni e trasformarle da soli, e usare il loro risolutore per trovare punti di intersezione.

Answer 3

È possibile utilizzare il metodo illustrato di seguito: https://math.stackexchange.com/questions/864070/how-to-determine-if-two-ellipse-have-at-least-one-intersection-point/864186#864186

In primo luogo si dovrebbe essere in grado di ridimensionare un'ellisse in una direzione. Per fare questo è necessario calcolare i coefficienti dell'ellisse come una sezione conica, ridimensionare e quindi recuperare i nuovi parametri geometrici dell'ellisse: centro, assi, angolo.

Quindi il problema si riduce a quello di trovare la distanza da un'ellisse all'origine. Per risolvere questo ultimo problema è necessario un po 'di iterazione. Ecco una possibile implementazione autonomo ...

from math import * 

def bisect(f,t_0,t_1,err=0.0001,max_iter=100): 
    iter = 0 
    ft_0 = f(t_0) 
    ft_1 = f(t_1) 
    assert ft_0*ft_1 <= 0.0 
    while True: 
     t = 0.5*(t_0+t_1) 
     ft = f(t) 
     if iter>=max_iter or ft<err: 
      return t 
     if ft * ft_0 <= 0.0: 
      t_1 = t 
      ft_1 = ft 
     else: 
      t_0 = t 
      ft_0 = ft 
     iter += 1 

class Ellipse(object): 
    def __init__(self,center,radius,angle=0.0): 
     assert len(center) == 2 
     assert len(radius) == 2 
     self.center = center 
     self.radius = radius 
     self.angle = angle 

    def distance_from_origin(self): 
     """                   
     Ellipse equation:                
     (x-center_x)^2/radius_x^2 + (y-center_y)^2/radius_y^2 = 1      
     x = center_x + radius_x * cos(t)            
     y = center_y + radius_y * sin(t)            
     """ 
     center = self.center 
     radius = self.radius 

     # rotate ellipse of -angle to become axis aligned        

     c,s = cos(self.angle),sin(self.angle) 
     center = (c * center[0] + s * center[1], 
        -s* center[0] + c * center[1]) 

     f = lambda t: (radius[1]*(center[1] + radius[1]*sin(t))*cos(t) - 
         radius[0]*(center[0] + radius[0]*cos(t))*sin(t)) 

     if center[0] > 0.0: 
      if center[1] > 0.0: 
       t_0, t_1 = -pi, -pi/2 
      else: 
       t_0, t_1 = pi/2, pi 
     else: 
      if center[1] > 0.0: 
       t_0, t_1 = -pi/2, 0 
      else: 
       t_0, t_1 = 0, pi/2 

     t = bisect(f,t_0,t_1) 
     x = center[0] + radius[0]*cos(t) 
     y = center[1] + radius[1]*sin(t) 
     return sqrt(x**2 + y**2) 

print Ellipse((1.0,-1.0),(2.0,0.5)).distance_from_origin()