essai Ray-/segment intersection pour le parallélisme et colinéarité ne BCZ de précision float en python

Question

Je suis en train de mettre en œuvre une fonction pour trouver des intersections Ray-/segment en python suivant les instructions de grandes Gareth Rees: https://stackoverflow.com/a/14318254/7235455 et https://stackoverflow.com/a/565282/7235455

Voilà ma fonction:

from math import radians, sin, cos 
import numpy as np 

def find_intersection(point0, theta, point1, point2): 
    # convert arguments to arrays: 
    p = np.array(point0, dtype=np.float) # ray origin 
    q = np.array(point1, dtype=np.float) # segment point 1 
    q2 = np.array(point2, dtype=np.float) # segment point 2 
    r = np.array((cos(theta),sin(theta))) # theta as vector (= ray as vector) 

    s = q2 - q # vector from point1 to point2 
    rxs = np.cross(r,s) 
    qpxs = np.cross(q-p,s) 
    qpxr = np.cross(q-p,r) 
    t = qpxs/rxs 
    u = qpxr/rxs 

    if rxs == 0 and qpxr == 0: 
     t0 = np.dot(q-p,r)/np.dot(r,r) 
     t1 = np.dot(t0+s,r)/np.dot(r,r) 
     return "collinear" 
    elif rxs == 0 and qpxr != 0: 
     return "parallel" 
    elif rxs != 0 and 0 <= t and 0 <= u and u <= 1: # removed t <= 1 since ray is inifinte 
     intersection = p+t*r 
     return "intersection is {0}".format(intersection) 
    else: 
     return None 

La fonction fonctionne très bien quand il y a une intersection. Mais il ne reconnaît pas le parallélisme ou la colinéarité, car les conditions rxs == 0 et qpxr == 0 ne sont pas (jamais?) Remplies. Exécuter par exemple:

p0 = (0.0,0.0) 
theta = radians(45.0) 
p1 = (1.0,1.0) 
p2 = (3.0,3.0) 

c = find_intersection(p0,theta,p1,p2) 

qui renvoie None. Ajout d'une déclaration d'impression pour RXS et qpxr avant que le cas du bloc donne

rxs = 2.22044604925e-16 qpxr = -1.11022302463e-16 

Ma conclusion est, la fonction ne parvient pas à attraper les conditions de la première instruction if en raison de problèmes de virgule flottante. 2.22044604925e-16 et -1.11022302463e-16 sont assez petits, mais malheureusement pas exactement 0. Je comprends que les flottants ne peuvent pas avoir une représentation exacte en binaire.

est ma conclusion correcte ou ai-je manqué quelque chose? Existe-t-il des idées de mise en œuvre pour éviter ce problème? Merci beaucoup!

Answer 1

Oui, votre conclusion est juste, problème de la stabilité numérique du prédicat « parallèle ».

Vous pouvez comparer les résultats avec un nombre faible (par exemple, eps=1.0E-9). Sa magnitude peut dépendre de la plage de coordonnées (notez que le produit croisé donne une zone de triangle doublée, donc la normalisation eps par MaxVecLen**2 semble raisonnable).

option Plus complexe, mais plus exactement - en utilisant des prédicats géométriques robustes comme these ones. Peut-être que les bibliothèques Python/NumPy pour la géométrie computationnelle contiennent des implémentations pour de telles opérations.

Answer 2

Il existe un moyen simple et sûr de résoudre ce problème.

Écrivez l'équation implicite du rayon (S(X, Y) = a X + b Y + c = 0). Lorsque vous branchez les coordonnées des extrémités du segment dans la fonction S, vous obtenez deux valeurs, soit S0 et S1. Si elles sont de signe opposé, il y a une intersection entre la ligne de support du rayon et le segment.

Dans ce cas, la position de l'intersection le long du segment est donnée par la valeur du paramètre, ce qui équivaut à

- S0/(S1 - S0). 

Cette expression jouit de la propriété d'être toujours calculable (condition qu'il y ait un changement de signe) et dans la plage [0, 1]. Il permet de calculer en toute sécurité le point d'intersection.

Pour ne sélectionner que les intersections qui se trouvent sur la demi-droite désirée (rayon), il suffit de calculer le signe S(Xo, Yo), à l'origine du rayon.

Cette procédure ne détecte pas les rayons parallèles ni colinéaires, mais il n'a pas d'importance. En tout cas, il produit des résultats sonores.