¿Cuál es la forma más eficiente de encontrar coordenadas barcéntricas?

45

En mi perfilador, encontrar coordenadas barcéntricas es aparentemente un cuello de botella. Estoy buscando hacerlo más eficiente.

Sigue el método en shirley , donde calculas el área de los triángulos formados incrustando el punto P dentro del triángulo.

bary

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, ¡pero estoy buscando uno más eficiente!

bobobobo
fuente
2
Tenga en cuenta que las soluciones más eficientes pueden ser las menos precisas.
Peter Taylor
Le sugiero que haga una prueba unitaria para llamar a este método ~ 100k veces (o algo similar) y mida el rendimiento. Puede escribir una prueba que garantice que sea menor que algún valor (p. Ej., 10s), o puede usarla simplemente para comparar la implementación antigua con la nueva.
cenizas999

Respuestas:

54

Transcrito de Detección de colisión en tiempo real de Christer Ericson (que, por cierto, es un libro excelente):

// 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 es efectivamente la regla de Cramer para resolver un sistema lineal. No será mucho más eficiente que esto: si esto sigue siendo un cuello de botella (y podría serlo: no parece que sea muy diferente en cuanto a cómputo que su algoritmo actual), probablemente necesitará encontrar otro lugar para ganar una aceleración

Tenga en cuenta que un número decente de valores aquí es independiente de p; se pueden almacenar en caché con el triángulo si es necesario.

John Calsbeek
fuente
77
El número de operaciones puede ser un arenque rojo. La forma en que dependen y los horarios es muy importante en las CPU modernas. siempre pruebe los supuestos y las "mejoras" de rendimiento.
Sean Middleditch
1
Las dos versiones en cuestión tienen una latencia casi idéntica en la ruta crítica, si solo está mirando operaciones matemáticas escalares. Lo que me gusta de este es que al pagar espacio por solo dos flotadores, puede eliminar una resta y una división de la ruta crítica. ¿ Eso vale la pena? Sólo una prueba de rendimiento sabe a ciencia cierta ...
John Calsbeek
1
Describe cómo consiguió esto en la página 137-138 con la sección "punto más cercano del triángulo al punto"
bobobobo
1
Nota menor: no hay argumento ppara esta función.
Bart
2
Nota de implementación menor: si los 3 puntos están uno encima del otro, obtendrá un error "dividir por 0", así que asegúrese de verificar ese caso en el código real.
frodo2975
9

La regla de Cramer debería ser la mejor manera de resolverlo. No soy un tipo gráfico, pero me preguntaba por qué en el libro Real-Time Collision Detection no hacen lo siguiente más simple:

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

Esto resuelve directamente el sistema lineal 2x2

v v0 + w v1 = v2

mientras que el método del libro resuelve el sistema

(v v0 + w v1) dot v0 = v2 dot v0
(v v0 + w v1) dot v1 = v2 dot v1
usuario5302
fuente
¿Su solución propuesta no hace suposiciones sobre la tercera .zdimensión () (específicamente, que no existe)?
Cornstalks
1
Este es el mejor método aquí si uno está trabajando en 2D. Solo una mejora menor: uno debe calcular el recíproco del denominador para usar dos multiplicaciones y una división en lugar de dos divisiones.
rubik
8

Ligeramente más rápido: calcula previamente el denominador y multiplica en lugar de dividir. Las divisiones son mucho más caras que las multiplicaciones.

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

En mi implementación, sin embargo, almacené en caché todas las variables independientes. Precalculo lo siguiente en el constructor:

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

Entonces el código final se ve así:

// 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
fuente
2

Usaría la solución que publicó John, pero usaría el SSS 4.2 dot intrínseco y sse rcpss intrínseco para la división, suponiendo que esté bien restringiéndose a Nehalem y procesos más nuevos y precisión limitada.

Alternativamente, podría calcular varias coordenadas barcéntricas a la vez utilizando sse o avx para una aceleración de 4 u 8x.

Crowley9
fuente
1

Puede convertir su problema 3D en un problema 2D proyectando uno de los planos alineados a los ejes y utilizando el método propuesto por user5302. Esto dará como resultado exactamente las mismas coordenadas barcéntricas siempre y cuando se asegure de que su triángulo no se proyecte en una línea. Lo mejor es proyectar en el plano alineado al eje que esté lo más cerca posible de la orientación de su triángulo. Esto evita problemas de co-linealidad y garantiza la máxima precisión.

En segundo lugar, puede calcular previamente el denominador y almacenarlo para cada triángulo. Esto guarda los cálculos posteriores.

Gert
fuente
1

Intenté copiar el código de @ NielW a C ++, pero no obtuve los resultados correctos.

Fue más fácil leer https://en.wikipedia.org/wiki/Barycentric_coordinate_system#Barycentric_coordinates_on_triangles y calcular el lambda1 / 2/3 como se indica allí (no se necesitan funciones vectoriales).

Si p (0..2) son los puntos del triángulo con x / y / z:

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

entonces las lambdas para un punto "punto" son

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;
usuario1712200
fuente
0

Para un punto dado N dentro del triángulo ABC, puede obtener el peso barcéntrico del punto C dividiendo el área del subtriángulo ABN por el área total del triángulo AB C.

Gandul
fuente