Existe-t-il un moyen itératif de calculer les rayons le long d’une ligne de balayage?

https://stackoverflow.com/questions/1423918

07-07-2019
|

Question

Je traite une série de points qui ont tous la même valeur Y, mais des valeurs X différentes. Je passe les points en incrémentant X de un. Par exemple, je pourrais avoir Y = 50 et X étant les entiers compris entre -30 et 30. Une partie de mon algorithme consiste à rechercher la distance à l'origine à partir de chaque point, puis à effectuer un traitement supplémentaire.

Après avoir établi le profil, j’ai constaté que l’appel sqrt dans le calcul de la distance prend beaucoup de temps. Existe-t-il un moyen itératif de calculer la distance?

En d'autres termes:

Je souhaite calculer efficacement: r[n] = sqrt(x[n]*x[n] + y*y)). Je peux sauvegarder des informations de l'itération précédente. Chaque itération change en incrémentant x, donc x[n] = x[n-1] + 1. Je ne peux pas utiliser les fonctions sqrt ou trig car elles sont trop lentes sauf au début de chaque ligne de balayage.

Je peux utiliser des approximations dans la mesure où elles sont suffisamment bonnes (moins de 0,1% d'erreur) et que les erreurs introduites sont lisses (je ne peux pas utiliser une table d'approximations précalculée).

Informations complémentaires: x et y sont toujours des entiers compris entre -150 et 150

Je vais essayer quelques idées demain et marquer la meilleure réponse en fonction de la plus rapide.

Résultats

J'ai fait quelques timings

Formule de distance: 16 ms / itération
La solution d'interpolation de Pete: 8 ms / itération
solution de pré-calcul wrang-wrang: 8 ms / itération

J'espérais que le test déciderait entre les deux, parce que j'aime les deux réponses. Je vais y aller avec Pete parce qu'il utilise moins de mémoire.

La solution

Juste pour avoir une idée, pour votre gamme y = 50, x = 0 donne r = 50 et y = 50, x = +/- 30 donne r ~ = 58.3. Vous voulez une approximation bonne pour +/- 0,1%, ou +/- 0,05 absolu. C'est beaucoup moins précis que la plupart des bibliothèques.

Deux approches approximatives - vous calculez r en vous basant sur l’interpolation de la valeur précédente ou utilisez quelques termes d’une série appropriée.

Interpolant à partir de r précédent

r = (x ² + y ²) ^1/2

dr / dx = 1/2. 2x. (x ² + y ²) ^-1/2 = x / r

    double r = 50;

    for ( int x = 0; x <= 30; ++x ) {

        double r_true = Math.sqrt ( 50*50 + x*x );

        System.out.printf ( "x: %d r_true: %f r_approx: %f error: %f%%\n", x, r, r_true, 100 * Math.abs ( r_true - r ) / r );

        r = r + ( x + 0.5 ) / r; 
    }

donne:

x: 0 r_true: 50.000000 r_approx: 50.000000 error: 0.000000%
x: 1 r_true: 50.010000 r_approx: 50.009999 error: 0.000002%
....
x: 29 r_true: 57.825065 r_approx: 57.801384 error: 0.040953%
x: 30 r_true: 58.335225 r_approx: 58.309519 error: 0.044065%

qui semble répondre à l'exigence d'erreur de 0,1%, je ne me suis donc pas soucié de coder la suivante, car cela demanderait un peu plus d'étapes de calcul.

Série tronquée

La série taylor pour sqrt (1 + x) pour x proche de zéro est

sqrt (1 + x) = 1 + 1/2 x - 1/8 x ² ... + (- 1/2) ^{n + 1} x ⁿ

En utilisant r = y sqrt (1 + (x / y) ²), alors vous recherchez un terme t = (- 1/2) ^{n + 1} 0.36 ⁿ avec une magnitude inférieure à 0.001, log (0.002) & Gt; n log (0,18) ou n > 3.6, donc prendre les termes à x ^ 4 devrait être Ok.

Autres conseils

Y=10000
Y2=Y*Y
for x=0..Y2 do
  D[x]=sqrt(Y2+x*x)

norm(x,y)=
  if (y==0) x
  else if (x>y) norm(y,x) 
  else {
     s=Y/y
     D[round(x*s)]/s
  }

Si vos coordonnées sont lisses, l’idée peut être étendue avec une interpolation linéaire. Pour plus de précision, augmentez Y.

L'idée est que s * (x, y) est sur la ligne y = Y, pour laquelle vous avez pré-calculé les distances. Obtenez la distance, puis divisez-la par S.

Je suppose que vous avez vraiment besoin de la distance et non de son carré.

Vous pourrez peut-être également trouver une implémentation générale du sqrt qui sacrifie un peu de précision pour la rapidité, mais j’ai du mal à imaginer que battre ce que la FPU peut faire.

Par interpolation linéaire, je souhaite remplacer D[round(x)] par:

f=floor(x)
a=x-f
D[f]*(1-a)+D[f+1]*a

Cela ne répond pas vraiment à votre question, mais peut aider ...

Les premières questions que je voudrais poser seraient:

& "ai-je besoin du sqrt du tout? &";
& "; Sinon, comment puis-je réduire le nombre de sqrts? &";
puis le vôtre: & "Puis-je remplacer les carrés restants par un calcul intelligent? &";

Je commencerais donc par:

Avez-vous besoin du rayon exact ou un rayon carré serait-il acceptable? Il y a approximatiosn rapide à sqrt, mais probablement pas assez précis pour votre spec.
Pouvez-vous traiter l'image à l'aide de quadrants ou de huitièmes en miroir? En traitant tous les pixels à la même valeur de rayon dans un lot, vous pouvez réduire le nombre de calculs de 8x.
Pouvez-vous précalculer les valeurs de rayon? Vous n'avez besoin que d'une table représentant un quart (ou éventuellement un huitième) de la taille de l'image que vous traitez, et il suffirait de la précalculer une fois, puis de la réutiliser pour de nombreuses exécutions de l'algorithme.

Donc, une mathématique intelligente peut ne pas être la solution la plus rapide.

Eh bien, il y a toujours des efforts pour optimiser votre place, le plus rapide que j'ai vu est le vieux tremblement automobile 3 sqrt:

http://betterexplained.com/articles/understanding- quakes-fast-inverse-square-root /

Cela dit, comme sqrt est non linéaire, vous ne pourrez pas effectuer d'interpolation linéaire simple le long de votre ligne pour obtenir votre résultat. La meilleure idée est d'utiliser une table de consultation, car cela vous donnera un accès extrêmement rapide aux données. Et, comme vous semblez itérer par des entiers entiers, une recherche dans un tableau devrait être extrêmement précise.

Eh bien, vous pouvez commencer par refléter autour de x = 0 (il suffit de calculer n > = 0 et de dupliquer ces résultats avec n < 0 correspondant). Après cela, je regarderais comment utiliser la dérivée sur sqrt (a ^ 2 + b ^ 2) (ou le péché correspondant) pour tirer parti de la constante dx.

Si cela n’est pas assez précis, puis-je vous dire que c’est un très bon travail pour SIMD, qui vous fournira une opération de racine carrée réciproque sur SSE et VMX (et le modèle de shader 2).

Ceci est en quelque sorte lié à un élément HAKMEM :

POINT 149 (Minsky): ALGORITHME DE CERCLE Voici une manière élégante de dessiner presque cercles sur un affichage en points:
NEW X = OLD X - epsilon * OLD Y
NEW Y = OLD Y + epsilon * NEW(!) X
Cela donne une ellipse très ronde   centré à l'origine avec sa taille   déterminé par le point initial.   epsilon détermine l’angle   la vitesse du point de circulation, et   affecte légèrement l'excentricité. Si   epsilon est une puissance de 2, alors nous ne le faisons pas   même besoin de multiplication, encore moins   racines carrées, sinus et cosinus! le   " cercle " sera parfaitement stable   parce que les points deviennent bientôt   périodique.

L’algorithme de cercle a été inventé par   erreur quand j'ai essayé de sauver un   inscrivez-vous dans un hack d'affichage! Ben Gurley   avait un hack d'affichage incroyable en utilisant seulement   environ six ou sept instructions, et   c'était une grande merveille. Mais c'était   fondamentalement orienté ligne. C'est arrivé   pour moi qu'il serait excitant de   avoir des courbes, et j'essayais d'obtenir un   affichage de courbe hack avec un minimum   instructions.

Licencié sous: CC-BY-SA avec attribution

Non affilié à StackOverflow