Qual é a maneira mais eficiente de encontrar coordenadas baricêntricas?

No meu perfilador, encontrar coordenadas barentêntricas é aparentemente um gargalo. Eu estou olhando para torná-lo mais eficiente.

Segue o método em shirley , onde você calcula a área dos triângulos formados incorporando o ponto P dentro do triângulo.

Código:

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

Este método funciona, mas estou procurando um mais eficiente!

Transcrito da Detecção de colisão em tempo real de Christer Ericson (que, aliás, é um excelente livro):

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

Esta é efetivamente a regra de Cramer para resolver um sistema linear. Você não será muito mais eficiente do que isso - se isso ainda for um gargalo (e pode ser: não parece muito diferente em termos de computação que o algoritmo atual), você provavelmente precisará encontrar outro lugar para ganhar uma aceleração.

Observe que um número decente de valores aqui é independente de p - eles podem ser armazenados em cache com o triângulo, se necessário.

A regra do Cramer deve ser a melhor maneira de resolvê-la. Eu não sou um cara gráfico, mas queria saber por que, no livro Detecção de colisão em tempo real, eles não fazem a seguinte coisa mais simples:

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

Isso resolve diretamente o sistema linear 2x2

v v0 + w v1 = v2

enquanto o método do livro resolve o sistema

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

Um pouco mais rápido: pré-calcule o denominador e multiplique em vez de dividir. As divisões são muito mais caras que as multiplicações.

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

Na minha implementação, no entanto, eu armazenei em cache todas as variáveis ​​independentes. Eu pré-calculei o seguinte no construtor:

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

Portanto, o código final fica assim:

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

Eu usaria a solução que John postou, mas usaria o SSS 4.2 ponto intrínseco e o sse rcpss intrínseco para a divisão, supondo que você esteja bem restringindo-se a Nehalem e processos mais recentes e precisão limitada.

Como alternativa, você pode calcular várias coordenadas barricêntricas de uma só vez usando sse ou avx para uma aceleração de 4 ou 8x.

Você pode converter seu problema 3D em um problema 2D projetando um dos planos alinhados ao eixo e usar o método proposto pelo usuário5302. Isso resultará exatamente nas mesmas coordenadas barricêntricas, desde que você verifique se o triângulo não se projeta em uma linha. O melhor é projetar no plano alinhado ao eixo o mais próximo possível da orientação da triagle. Isso evita problemas de co-linearidade e garante a máxima precisão.

Em segundo lugar, você pode pré-calcular o denominador e armazená-lo para cada triângulo. Isso salva os cálculos posteriormente.

Tentei copiar o código do @ NielW para C ++, mas não obtive resultados corretos.

Era mais fácil ler https://en.wikipedia.org/wiki/Barycentric_coordinate_system#Barycentric_coordinates_on_triangles e calcular o lambda1 / 2/3 conforme indicado lá (não são necessárias funções vetoriais).

Se p (0..2) são os pontos do triângulo com x / y / z:

Pré-cálculo para triângulo:

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

então as lambdas para um Point "point" são

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;

Para um dado ponto N dentro do triângulo ABC, você pode obter o peso barcentric do ponto C dividindo a área do sub-triângulo ABN pela área total do triângulo AB C.