What's the most efficient way to find barycentric coordinates?

Question

In my profiler, finding barycentric coordinates is apparently somewhat of a bottleneck. I am looking to make it more efficient.

It follows the method in shirley, where you compute the area of the triangles formed by embedding the point P inside the triangle.

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

This method works, but I'm looking for a more efficient one!

Answer 1

Transcribed from Christer Ericson's Real-Time Collision Detection (which, incidentally, is an excellent book):

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

This is effectively Cramer's rule for solving a linear system. You will not get much more efficient than this—if this is still a bottleneck (and it might be: it doesn't look like it's much different computation-wise than your current algorithm), you'll probably need to find some other place to gain a speedup.

Note that a decent number of values here are independent of p—they can be cached with the triangle if necessary.

Answer 2

The Cramer's rule should be the best way to solve it. I am not a graphic guy, but I was wondering why in the book Real-Time Collision Detection they doesn't do the following simpler thing:

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

This directly solves the 2x2 linear system

v v0 + w v1 = v2

while the method from the book solves the system

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

Answer 3

Slightly faster: Precompute the denominator, and multiply instead of divide. Divisions are much more expensive than multiplications.

// 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 my implementation, however, I cached all of the independent variables. I pre-calc the following in the constructor:

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

So the final code looks like this:

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

Answer 4

I tried to copy @NielW's code to C++, but I didn't get correct results.

It was easier to read https://en.wikipedia.org/wiki/Barycentric_coordinate_system#Barycentric_coordinates_on_triangles and calculate the lambda1/2/3 as given there (no vector functions needed).

If p(0..2) are the Points of the triangle with x/y/z:

Precalc for triangle:

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

then the lambdas for a Point "point" are

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;

Answer 5

I would use the solution that John posted, but I would use the SSS 4.2 dot intrinsic and sse rcpss intrinsic forthe divide, assuming you are ok restricting yourself to Nehalem and newer processes and limited precision.

Alternatively you could compute several barycentric coordinates at once using sse or avx for a 4 or 8x speedup.

Answer 6

You can convert your 3D problem into a 2D problem by projecting one of the axis-aligned planes and use the method proposed by user5302. This will result in exactly the same barycentric coordinates as long as you make sure your triangle does not project into a line. Best is to project to the axis-aligned plane that is as close as possible to the orientation of your triagle. This avoid co-linearity problems and ensure maximum accuracy.

Secondly you can pre-compute the denominator and store it for each triangle. This saves computations afterwards.

Answer 7

For a given point N inside triangle A B C, you can get the barycentric weight of point C by dividing the area of subtriangle A B N by the total area of triangle A B C.

Answer 8

It may not be the fastest answer, but it's definitely one of the cleanest: You can convert to homogenous coordinates in order to use Cramer's rule cleanly for each barycentric coordinate. Please ignore the strange syntax, it is a custom preprocessor I wrote over OpenCL C.

VEC3 barycentric(VEC3 P, Triangle3f tri){

    VEC3 A = VEC3::diff(_VEC3(tri.A),P);
    VEC3 B = VEC3::diff(_VEC3(tri.B),P);
    VEC3 C = VEC3::diff(_VEC3(tri.C),P);
    
    float Ax = A.unwrap()..x;
    float Ay = A.unwrap()..y;
    float Az = A.unwrap()..z;
    
    float Bx = B.unwrap()..x;
    float By = B.unwrap()..y;
    float Bz = B.unwrap()..z;
    
    float Cx = C.unwrap()..x;
    float Cy = C.unwrap()..y;
    float Cz = C.unwrap()..z;
    
    float Px = 0.0;
    float Py = 0.0;
    float Pz = 0.0;
    
    MAT4 M = __MAT4(
    
        Ax,Bx,Cx,1,
        Ay,By,Cy,1,
        Az,Bz,Cz,1,
        1, 1, 1, 0
    
    );
    
    float d = MAT4::det(M).scalarUnwrap();
    
    if(fabs(d)<DET_EPSILON){
        return VEC3::zero();
    }
    
    MAT4 M1 = __MAT4(
    
        Px,Bx,Cx,1,
        Py,By,Cy,1,
        Pz,Bz,Cz,1,
        1, 1, 1, 0
    
    );
    
    MAT4 M2 = __MAT4(
    
        Ax,Px,Cx,1,
        Ay,Py,Cy,1,
        Az,Pz,Cz,1,
        1, 1, 1, 0
    
    );
    
    MAT4 M3 = __MAT4(
    
        Ax,Bx,Px,1,
        Ay,By,Py,1,
        Az,Bz,Pz,1,
        1, 1, 1, 0
    
    );
    
    return VEC3::fromcoords(
        MAT4::det(M1).scalarUnwrap() / d,
        MAT4::det(M2).scalarUnwrap() / d,
        MAT4::det(M3).scalarUnwrap() / d
    );

}