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

45

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.

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

Diese Methode funktioniert, aber ich suche eine effizientere!

Bobobobo
quelle
2
Beachten Sie, dass die effizientesten Lösungen möglicherweise die ungenauesten sind.
Peter Taylor
Ich schlage vor, dass Sie einen Komponententest durchführen, um diese Methode ~ 100.000 Mal (oder etwas Ähnliches) aufzurufen und die Leistung zu messen. Sie können einen Test schreiben, der sicherstellt, dass er unter einem bestimmten Wert liegt (z. B. 10 Sekunden), oder Sie können ihn einfach verwenden, um alte oder neue Implementierungen zu vergleichen.
Asche999

Antworten:

54

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.

John Calsbeek
quelle
7
Anzahl der Operationen kann ein roter Hering sein. Wie abhängig und zeitlich abhängig sie sind, hängt stark von modernen CPUs ab. Testen Sie immer Annahmen und Leistungsverbesserungen.
Sean Middleditch
1
Die beiden fraglichen Versionen weisen auf dem kritischen Pfad eine nahezu identische Latenz auf, wenn Sie nur skalare mathematische Operationen betrachten. Das, was ich an diesem mag, ist, dass man durch das Bezahlen von Platz für nur zwei Schwimmer eine Subtraktion und eine Division vom kritischen Pfad entfernen kann. Lohnt sich das ? Nur ein Leistungstest weiß es mit Sicherheit ...
John Calsbeek
1
Wie er dazu kam, beschreibt er auf Seite 137-138 mit dem Abschnitt "Punkt am nächsten zum Punkt"
Bobobobo
1
Kleiner Hinweis: Es gibt kein Argument pfür diese Funktion.
Bart
2
Kleiner Hinweis zur Implementierung: Wenn alle 3 Punkte übereinander liegen, wird der Fehler "Teilen durch 0" angezeigt. Überprüfen Sie daher den tatsächlichen Code auf diesen Fall.
Frodo2975
9

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
user5302
quelle
Geht Ihre vorgeschlagene Lösung nicht von der dritten .zDimension () aus (insbesondere davon, dass sie nicht existiert)?
Cornstalks
1
Dies ist die beste Methode, wenn Sie in 2D arbeiten. Nur eine kleine Verbesserung: Man sollte den Kehrwert des Nenners berechnen, um zwei Multiplikationen und eine Division anstelle von zwei Divisionen zu verwenden.
Rubik
8

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;
}
NielW
quelle
2

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.

Crowley9
quelle
1

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.

Gert
quelle
1

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;
user1712200
quelle
0

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.

Dodger
quelle