Jaki jest najskuteczniejszy sposób znalezienia współrzędnych barycentrycznych?

W moim profilu znajdowanie współrzędnych barycentrycznych jest w pewnym sensie wąskim gardłem. Chcę sprawić, by był bardziej wydajny.

Podąża za metodą Shirleya , w której oblicza się powierzchnię trójkątów utworzonych przez osadzenie punktu P wewnątrz trójkąta.

Kod:

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

Ta metoda działa, ale szukam bardziej wydajnej!

Przepisywane z Christer Ericson's Real-Time Collision Detection (która, nawiasem mówiąc, jest doskonałą książką):

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

Jest to skutecznie reguła Cramera polegająca na rozwiązaniu układu liniowego. Nie zyskasz o wiele bardziej wydajnego niż to - jeśli nadal jest to wąskie gardło (a może być: nie wygląda na to, że różni się znacznie obliczeniami od obecnego algorytmu), prawdopodobnie będziesz musiał znaleźć inne miejsce aby przyspieszyć.

Zauważ, że przyzwoita liczba wartości jest tutaj niezależna od p - w razie potrzeby można je buforować za pomocą trójkąta.

Reguła Cramera powinna być najlepszym sposobem na jego rozwiązanie. Nie jestem grafikiem, ale zastanawiałem się, dlaczego w książce Wykrywanie kolizji w czasie rzeczywistym nie robią następującej prostszej rzeczy:

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

To bezpośrednio rozwiązuje układ liniowy 2x2

v v0 + w v1 = v2

podczas gdy metoda z książki rozwiązuje system

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

Nieco szybciej: oblicz wstępnie mianownik i pomnóż zamiast dzielić. Podziały są znacznie droższe niż mnożenie.

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

Jednak w mojej implementacji buforowałem wszystkie zmienne niezależne. Obliczam wstępnie w konstruktorze:

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

Ostateczny kod wygląda następująco:

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

Użyłbym rozwiązania opublikowanego przez Johna, ale użyłbym wewnętrznego SSS 4.2 i wewnętrznego rcpss dla podziału, zakładając, że nie przeszkadzasz sobie w Nehalem i nowszych procesach i ograniczonej precyzji.

Alternatywnie możesz obliczyć kilka współrzędnych barycentrycznych jednocześnie za pomocą sse lub avx dla przyspieszenia 4 lub 8x.

Możesz przekształcić swój problem 3D w problem 2D, rzutując jedną z płaszczyzn wyrównanych do osi i używając metody zaproponowanej przez user5302. Spowoduje to uzyskanie dokładnie takich samych współrzędnych barycentrycznych, o ile upewnisz się, że trójkąt nie wystaje w linię. Najlepiej jest rzutować na płaszczyznę wyrównaną do osi, która jest jak najbliżej orientacji twojego triagla. Pozwala to uniknąć problemów z współliniowością i zapewnia maksymalną dokładność.

Po drugie, możesz wstępnie obliczyć mianownik i zapisać go dla każdego trójkąta. Zapisuje to później obliczenia.

Próbowałem skopiować kod @ NielW do C ++, ale nie otrzymałem poprawnych wyników.

Łatwiej było przeczytać https://en.wikipedia.org/wiki/Barycentric_coordinate_system#Barycentric_coordinates_on_triangles i obliczyć lambda1 / 2/3 jak tam podano (żadne funkcje wektorowe nie są potrzebne).

Jeśli p (0..2) to punkty trójkąta o x / y / z:

Precalc dla trójkąta:

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

wtedy są lambdami dla „punktu” Punktu

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;

Dla danego punktu N wewnątrz trójkąta ABC można uzyskać masę barocentryczną punktu C, dzieląc obszar podtrójkąta ABN przez całkowitą powierzchnię trójkąta AB C.