Quel est le moyen le plus efficace de trouver des coordonnées barycentriques?

45

Dans mon profileur, trouver des coordonnées barycentriques est apparemment un peu un goulot d'étranglement. Je cherche à le rendre plus efficace.

Il suit la méthode de Shirley , où vous calculez l'aire des triangles formés en incorporant le point P à l'intérieur du triangle.

bary

Code:

Vector Triangle::getBarycentricCoordinatesAt( const Vector & P ) const
{
  Vector bary ;

  // The area of a triangle is 
  real areaABC = DOT( normal, CROSS( (b - a), (c - a) )  ) ;
  real areaPBC = DOT( normal, CROSS( (b - P), (c - P) )  ) ;
  real areaPCA = DOT( normal, CROSS( (c - P), (a - P) )  ) ;

  bary.x = areaPBC / areaABC ; // alpha
  bary.y = areaPCA / areaABC ; // beta
  bary.z = 1.0f - bary.x - bary.y ; // gamma

  return bary ;
}

Cette méthode fonctionne, mais je cherche une méthode plus efficace!

bobobobo
la source
2
Attention, les solutions les plus efficaces peuvent être les moins précises.
Peter Taylor
Je vous suggère de faire un test unitaire pour appeler cette méthode environ 100 000 fois (ou quelque chose de similaire) et mesurer les performances. Vous pouvez écrire un test qui garantit qu'elle est inférieure à une valeur (par exemple, 10s), ou vous pouvez l'utiliser simplement pour comparer une ancienne implémentation à une nouvelle.
ashes999

Réponses:

54

Transcrit à partir de la détection de collision en temps réel de Christer Ericson (qui est d'ailleurs un excellent livre):

// Compute barycentric coordinates (u, v, w) for
// point p with respect to triangle (a, b, c)
void Barycentric(Point p, Point a, Point b, Point c, float &u, float &v, float &w)
{
    Vector v0 = b - a, v1 = c - a, v2 = p - a;
    float d00 = Dot(v0, v0);
    float d01 = Dot(v0, v1);
    float d11 = Dot(v1, v1);
    float d20 = Dot(v2, v0);
    float d21 = Dot(v2, v1);
    float denom = d00 * d11 - d01 * d01;
    v = (d11 * d20 - d01 * d21) / denom;
    w = (d00 * d21 - d01 * d20) / denom;
    u = 1.0f - v - w;
}

C'est effectivement la règle de Cramer pour résoudre un système linéaire. Vous ne serez pas beaucoup plus efficace que cela - si cela reste un goulot d'étranglement (et cela pourrait être le cas: cela ne semble pas être très différent du point de vue du calcul par rapport à votre algorithme actuel), vous devrez probablement trouver un autre emplacement. pour gagner du temps.

Notez qu'un nombre décent de valeurs ici sont indépendantes de p - elles peuvent être mises en cache avec le triangle si nécessaire.

John Calsbeek
la source
7
Le nombre d'opérations peut être un hareng rouge. Leur dépendance et leur horaire comptent beaucoup pour les processeurs modernes. testez toujours les hypothèses et les «améliorations» des performances.
Sean Middleditch
1
Les deux versions en question ont une latence presque identique sur le chemin critique, si vous ne regardez que des opérations mathématiques scalaires. Ce que j’aime dans ce cas-là, c’est qu’en payant de l’espace pour seulement deux flottants, vous pouvez supprimer une soustraction et une division du chemin critique. Est-ce que ça vaut le coup? Seul un test de performance en est certain…
John Calsbeek
1
Il décrit comment il l'a obtenu aux pages 137-138 avec la section "Point le plus proche d'un triangle à l'autre"
bobobobo
1
Note mineure: il n'y a pas d'argument ppour cette fonction.
Bart
2
Note mineure d'implémentation: Si les 3 points sont superposés, vous obtiendrez une erreur "division par 0", assurez-vous donc de vérifier ce cas dans le code réel.
vendredi
9

La règle de Cramer devrait être le meilleur moyen de le résoudre. Je ne suis pas un graphiste, mais je me demandais pourquoi dans le livre Détection de collision en temps réel, ils ne font pas la chose plus simple suivante:

// Compute barycentric coordinates (u, v, w) for
// point p with respect to triangle (a, b, c)
void Barycentric(Point p, Point a, Point b, Point c, float &u, float &v, float &w)
{
    Vector v0 = b - a, v1 = c - a, v2 = p - a;
    den = v0.x * v1.y - v1.x * v0.y;
    v = (v2.x * v1.y - v1.x * v2.y) / den;
    w = (v0.x * v2.y - v2.x * v0.y) / den;
    u = 1.0f - v - w;
}

Cela résout directement le système 2x2 linéaire

v v0 + w v1 = v2

alors que la méthode du livre résout le système

(v v0 + w v1) dot v0 = v2 dot v0
(v v0 + w v1) dot v1 = v2 dot v1
utilisateur5302
la source
Votre solution proposée ne suppose-t-elle pas la troisième .zdimension (en particulier, le fait qu’elle n’existe pas)?
Cornstalks
1
C'est la meilleure méthode ici si l'on travaille en 2D. Juste une amélioration mineure: il faut calculer l'inverse du dénominateur afin d'utiliser deux multiplications et une division au lieu de deux divisions.
rubik
8

Légèrement plus rapide: Précalculez le dénominateur et multipliez au lieu de diviser. Les divisions sont beaucoup plus chères que les multiplications.

// Compute barycentric coordinates (u, v, w) for
// point p with respect to triangle (a, b, c)
void Barycentric(Point a, Point b, Point c, float &u, float &v, float &w)
{
    Vector v0 = b - a, v1 = c - a, v2 = p - a;
    float d00 = Dot(v0, v0);
    float d01 = Dot(v0, v1);
    float d11 = Dot(v1, v1);
    float d20 = Dot(v2, v0);
    float d21 = Dot(v2, v1);
    float invDenom = 1.0 / (d00 * d11 - d01 * d01);
    v = (d11 * d20 - d01 * d21) * invDenom;
    w = (d00 * d21 - d01 * d20) * invDenom;
    u = 1.0f - v - w;
}

Dans mon implémentation, cependant, j'ai mis en cache toutes les variables indépendantes. J'ai précalculé ce qui suit dans le constructeur:

Vector v0;
Vector v1;
float d00;
float d01;
float d11;
float invDenom;

Donc, le code final ressemble à ceci:

// Compute barycentric coordinates (u, v, w) for
// point p with respect to triangle (a, b, c)
void Barycentric(Point a, Point b, Point c, float &u, float &v, float &w)
{
    Vector v2 = p - a;
    float d20 = Dot(v2, v0);
    float d21 = Dot(v2, v1);
    v = (d11 * d20 - d01 * d21) * invDenom;
    w = (d00 * d21 - d01 * d20) * invDenom;
    u = 1.0f - v - w;
}
NielW
la source
2

J'utiliserais la solution publiée par John, mais j'utiliserais les points intrinsèques SSS 4.2 et intrinsèques pour la division, en supposant que vous vous en tenez à vous limiter à Nehalem et aux processus plus récents et à une précision limitée.

Alternativement, vous pouvez calculer plusieurs coordonnées barycentriques à la fois en utilisant sse ou avx pour une accélération de 4 ou 8x.

Crowley9
la source
1

Vous pouvez convertir votre problème 3D en problème 2D en projetant l’un des plans alignés sur l’axe et en utilisant la méthode proposée par user5302. Cela donnera exactement les mêmes coordonnées barycentriques tant que vous vous assurez que votre triangle ne se projette pas dans une ligne. Le mieux est de projeter sur le plan aligné sur l’axe le plus proche possible de l’orientation de votre triagle. Cela évite les problèmes de colinéarité et garantit une précision maximale.

Deuxièmement, vous pouvez pré-calculer le dénominateur et le stocker pour chaque triangle. Cela enregistre les calculs par la suite.

Gert
la source
1

J'ai essayé de copier le code de @ NielW en C ++, mais je n'ai pas obtenu de résultats corrects.

Il était plus facile de lire https://en.wikipedia.org/wiki/Barycentric_coordinate_system#Barycentric_coordinates_on_triangles et de calculer le lambda1 / 2/3 comme indiqué (aucune fonction vectorielle requise).

Si p (0..2) sont les points du triangle avec x / y / z:

Précalc pour triangle:

double invDET = 1./((p(1).y-p(2).y) * (p(0).x-p(2).x) + 
                   (p(2).x-p(1).x) * (p(0).y-p(2).y));

alors les lambdas pour un point "point" sont

double l1 = ((p(1).y-p(2).y) * (point.x-p(2).x) + (p(2).x-p(1).x) * (point.y-p(2).y)) * invDET; 
double l2 = ((p(2).y-p(0).y) * (point.x-p(2).x) + (p(0).x-p(2).x) * (point.y-p(2).y)) * invDET; 
double l3 = 1. - l1 - l2;
utilisateur1712200
la source
0

Pour un point N donné dans le triangle ABC, vous pouvez obtenir le poids barycentrique du point C en divisant l'aire du sous-angle ABN par l'aire totale du triangle AB C.

Dodger
la source