Was ist der effizienteste Weg, um baryzentrische Koordinaten zu finden?

In meinem Profiler ist das Finden von Schwerpunktkoordinaten anscheinend ein gewisser Engpass. Ich versuche es effizienter zu machen.

Es folgt der Methode in Shirley , bei der Sie die Fläche der durch Einbetten des Punkts P in das Dreieck gebildeten Dreiecke berechnen.

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 ;
}

Diese Methode funktioniert, aber ich suche eine effizientere!

Transkribiert von Christer Ericsons Echtzeit-Kollisionserkennung (die übrigens ein ausgezeichnetes Buch ist):

// 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;
}

Dies ist effektiv die Cramer-Regel für die Lösung eines linearen Systems. Sie werden nicht viel effizienter als dies sein - wenn dies immer noch ein Engpass ist (und es könnte sein, dass es in Bezug auf die Berechnung nicht viel anders aussieht als Ihr aktueller Algorithmus), müssen Sie wahrscheinlich einen anderen Ort finden eine Beschleunigung erreichen.

Beachten Sie, dass hier eine anständige Anzahl von Werten unabhängig von p ist - sie können bei Bedarf mit dem Dreieck zwischengespeichert werden.

Die Cramer-Regel sollte der beste Weg sein, um es zu lösen. Ich bin kein Grafiker, aber ich habe mich gefragt, warum in dem Buch Real-Time Collision Detection nicht die folgende einfache Sache gemacht wird:

// 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;
}

Dies löst direkt das 2x2-Linearsystem

v v0 + w v1 = v2

während die Methode aus dem Buch das System löst

(v v0 + w v1) dot v0 = v2 dot v0
(v v0 + w v1) dot v1 = v2 dot v1

Etwas schneller: Den Nenner vorberechnen und multiplizieren statt dividieren. Divisionen sind viel teurer als Multiplikationen.

// 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;
}

In meiner Implementierung habe ich jedoch alle unabhängigen Variablen zwischengespeichert. Ich habe folgendes im Konstruktor vorberechnet:

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

Der endgültige Code sieht also so aus:

// 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;
}

Ich würde die von John veröffentlichte Lösung verwenden, aber ich würde für die Teilung die intrinsische Punktzahl SSS 4.2 und die intrinsische Punktzahl sse rcpss verwenden, vorausgesetzt, Sie beschränken sich in Ordnung auf Nehalem und neuere Prozesse und beschränkte Präzision.

Alternativ können Sie mit sse oder avx mehrere Schwerpunktkoordinaten gleichzeitig für eine 4- oder 8-fache Beschleunigung berechnen.

Sie können Ihr 3D-Problem in ein 2D-Problem umwandeln, indem Sie eine der achsenausgerichteten Ebenen projizieren und die von user5302 vorgeschlagene Methode anwenden. Dies führt zu genau den gleichen Schwerpunktkoordinaten, solange Sie sicherstellen, dass Ihr Dreieck nicht in eine Linie hineinragt. Am besten projizieren Sie auf die achsenausgerichtete Ebene, die so nah wie möglich an der Ausrichtung Ihres Triagles liegt. Dies vermeidet Co-Linearitätsprobleme und gewährleistet maximale Genauigkeit.

Zweitens können Sie den Nenner vorberechnen und für jedes Dreieck speichern. Dies spart später Berechnungen.

Ich habe versucht, den Code von @ NielW nach C ++ zu kopieren, aber ich habe keine korrekten Ergebnisse erhalten.

Es war einfacher, https://en.wikipedia.org/wiki/Barycentric_coordinate_system#Barycentric_coordinates_on_triangles zu lesen und das Lambda1 / 2/3 wie dort angegeben zu berechnen (keine Vektorfunktionen erforderlich).

Wenn p (0..2) die Punkte des Dreiecks mit x / y / z sind:

Precalc für Dreieck:

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));

dann sind die Lambdas für einen Punkt "Punkt"

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;

Für einen gegebenen Punkt N innerhalb des Dreiecks ABC können Sie das Schwerpunktgewicht von Punkt C erhalten, indem Sie die Fläche des Subtriangels ABN durch die Gesamtfläche des Dreiecks AB C dividieren.