Calcul d'ordonnée d'intersection de lignes le plus précis avec des flottants?

Question

Je calcule l'ordonnée y d'un point sur une ligne à une abscisse x donnée. La ligne est définie par ses deux coordonnées (x0, y0) (x1, y1). Les coordonnées des points d'extrémité sont des flottants et le calcul doit être effectué avec une précision flottante pour une utilisation dans GPU.

Les mathématiques, et donc l'implémentation naïve, sont triviales. Soit t = (x - x0)/(x1 - x0), alors y = (1 - t) * y0 + t * y1 = y0 + t * (y1 - y0).

Le problème est lorsque x1 - x0 est petit. Le résultat introduira une erreur d'annulation. Quand il est combiné avec celui de x - x0, dans la division j'attends une erreur significative dans t.

La question est de savoir s'il existe une autre façon de déterminer y avec une meilleure précision?

est-ce que je devrais calculer (x - x0) * (y1 - y0) en premier et diviser par (x1 - x0) après?

La différence y1 - y0 sera toujours grande.

Answer 1

Dans une large mesure, votre problème sous-jacent est fondamental. Quand (x1-x0) est petit, cela signifie qu'il y a seulement quelques bits dans la mantisse de x1 et x0 qui diffèrent. Et par extension, il n'y a qu'un nombre limité de flotteurs entre x0 et x1. Par exemple.si seulement les 4 bits inférieurs de la mantisse diffèrent, il y a au plus 14 valeurs entre eux.

Dans votre meilleur algorithme, le terme t représente ces bits inférieurs. Et pour continuer ou exemple, si x0 et x1 diffèrent de 4 bits, alors t ne peut prendre que 16 valeurs. Le calcul de ces valeurs possibles est assez robuste. Que vous calculiez 3E0/14E0 ou ​​3E-12/14E-12, le résultat sera proche de la valeur mathématique de 3/14.

Votre formule a l'avantage supplémentaire d'avoir y0 < = y < = y1, puisque 0 < = t < = 1

(je suppose que vous en savez assez sur les représentations de flotteur, et par conséquent « (x1 -x0) est petit "signifie vraiment" petit, par rapport aux valeurs de x1 et x0 eux-mêmes "Une différence de 1E-1 est petit quand x0 = 1E3 mais grand si x0 = 1E-6)

Answer 2

Si vous avez la possibilité de le faire, vous pouvez introduire deux cas dans votre calcul, en fonction de abs (x1-x0) < abs (y1-y0). Dans le cas vertical abs (x1-x0) < abs (y1-y0), calculer x à partir de y au lieu de y à partir de x.

EDIT. Une autre possibilité serait d'obtenir le résultat bit par bit en utilisant une variante de recherche dichotomique. Ce sera plus lent, mais peut améliorer le résultat dans des cas extrêmes.

// Input is X 
xmin = min(x0,x1); 
xmax = max(x0,x1); 
ymin = min(y0,y1); 
ymax = max(y0,y1); 
for (int i=0;i<20;i++) // get 20 bits in result 
{ 
    xmid = (xmin+xmax)*0.5; 
    ymid = (ymin+ymax)*0.5; 
    if (x < xmid) { xmax = xmid; ymax = ymid; } // first half 
    else { xmin = xmid; ymin = ymid; } // second half 
} 
// Output is some value in [ymin,ymax] 
Y = ymin; 

Answer 3

Si vos données sources sont déjà flottantes, vous avez déjà une imprécision fondamentale.

Pour expliquer plus en détail, imaginez si vous le faites graphiquement. Vous avez une feuille 2D de papier millimétré, et 2 points marqués.

Cas 1: Ces points sont très précis, et ont été marqués avec un crayon très pointu. Il est facile de dessiner la ligne qui les relie, et facile d'y obtenir x (ou vice versa). Cas 2: Ces pointes ont été marquées d'un gros feutre, comme un marqueur de bingo. Clairement la ligne que vous dessinez sera moins précise. Passez-vous par le centre des spots? Le bord supérieur? Le bord inférieur? En haut de l'un, en bas de l'autre? Clairement, il y a beaucoup d'options différentes. Si les deux points sont proches l'un de l'autre, la variation sera encore plus grande. Les flotteurs ont un certain niveau d'imprécision inhérent à eux, en raison de la façon dont ils représentent les nombres, dans la mesure où ils correspondent plus au cas 2 qu'au cas 1 (ce qui pourrait être l'équivalent d'un librray de précision arbitraire). Aucun algorithme dans le monde ne peut compenser cela. Données imprécises dans, Données imprécises sur

Answer 4

Vérifiez que la distance entre x0 et x1 est petite, c'est-à-dire que fabs (x1 - x0) < eps. Ensuite, la ligne est parallèle à l'axe y du système de coordonnées, c'est-à-dire que vous ne pouvez pas calculer les valeurs y de cette ligne en fonction de sur x. Vous avez infini de nombreuses valeurs y et vous devez donc traiter ce cas différemment.

Answer 5

quelque chose comme:

t = sign * power2 (sqrt (abs(x - x0))/ sqrt (abs(x1 - x0))) 

L'idée est d'utiliser une formule mathématique équivalente dans laquelle low (x1-x0) a moins d'effet. (Je ne sais pas si celui que j'ai écrit correspond à ce critère)

Answer 6

J'ai mis en œuvre un programme de référence pour comparer l'effet de l'expression différente.

J'ai calculé y en utilisant une double précision, puis j'ai calculé en utilisant une seule précision avec des expressions différentes.

Voici l'expression testée:

inline double getYDbl(double x, double x0, double y0, double x1, double y1) 
{ 
    double const t = (x - x0)/(x1 - x0); 
    return y0 + t*(y1 - y0); 
} 

inline float getYFlt1(float x, float x0, float y0, float x1, float y1) 
{ 
    double const t = (x - x0)/(x1 - x0); 
    return y0 + t*(y1 - y0); 
} 

inline float getYFlt2(float x, float x0, float y0, float x1, float y1) 
{ 
    double const t = (x - x0)*(y1 - y0); 
    return y0 + t/(x1 - x0); 
} 

inline float getYFlt3(float x, float x0, float y0, float x1, float y1) 
{ 
    double const t = (y1 - y0)/(x1 - x0); 
    return y0 + t*(x - x0); 
} 

inline float getYFlt4(float x, float x0, float y0, float x1, float y1) 
{ 
    double const t = (x1 - x0)/(y1 - y0); 
    return y0 + (x - x0)/t; 
} 

I calculé la moyenne et stdDev de la différence entre le résultat double précision et le résultat de simple précision.

Le résultat est qu'il n'y en a aucun en moyenne plus de 1000 et 10K ensembles de valeurs aléatoires. J'ai utilisé le compilateur icc avec et sans optimisation ainsi que g ++.

Notez que j'ai dû utiliser la fonction isnan() pour filtrer les valeurs fausses. Je soupçonne que ceux-ci résultent d'un débordement dans la différence ou la division.

Je ne sais pas si les compilateurs réorganisent l'expression. Quoi qu'il en soit, la conclusion de ce test est que les réarrangements de l'expression ci-dessus n'ont aucun effet sur la précision de calcul. L'erreur reste la même (en moyenne).

Answer 7

Vous pouvez jeter un coup d'œil aux sources "QLine" de Qt (si je me souviens bien); ils ont implémenté un algorithme de détermination des intersections tiré d'un des livres "Graphics Gems" (la référence doit être dans les commentaires du code, le livre était sur EDonkey il y a quelques années), qui, à son tour, a des garanties sur l'applicabilité résolution d'écran donnée lorsque les calculs sont effectués avec une largeur de bit donnée (ils utilisent l'arithmétique à virgule fixe si je ne me trompe pas).

Answer 8

Comme indiqué MSalters, le problème est déjà dans les données d'origine.

L'interpolation/extrapolation nécessite la pente, qui a déjà une faible précision dans les conditions données (pire pour les segments de ligne très courts éloignés de l'origine).

Choix de l'algorithme canot retrouver cette perte de précision.Mon intuition est que l'ordre d'évaluation différent ne changera pas les choses, car l'erreur est introduite par les soustractions, pas la devision.

---

Idée:
Si vous avez des données plus précises lorsque les lignes sont générées, vous pouvez changer la représentation de ((x0, y0), (x1, y1)) à (x0, y0, angle, longueur). Vous pouvez stocker l'angle ou la pente, la pente a un poteau, mais l'angle nécessite des fonctions trigonométriques ... moche.

Bien sûr, cela ne fonctionnera pas si vous avez besoin du point de terminaison fréquemment, et vous avez tellement de lignes que vous ne pouvez pas stocker de données supplémentaires, je n'en ai aucune idée. Mais peut-être existe-t-il une autre représentation qui fonctionne bien pour vos besoins. Les doubles ont suffisamment de résolution dans la plupart des situations, mais cela doublerait également le jeu de travail. Pouvez-vous utiliser des fonctions trigonométriques comme sin ou cos?