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

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!

• Beware that the most efficient solutions may be the least accurate. Feb 12, 2012 at 13:51
• I suggest you make a unit test to call this method ~100k times (or something similar) and measure the performance. You can write a test that ensures it's less than some value (eg. 10s), or you can use it simply to benchmark old vs. new implementation. Feb 13, 2012 at 3:08

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.

• # of operations can be a red herring. How they're dependent and schedules matters a lot on modern CPUs. always test assumptions and performance "improvements." Feb 14, 2013 at 22:16
• The two versions in question have almost identical latency on the critical path, if you're only looking at scalar math ops. The thing I like about this one is that by paying space for merely two floats, you can shave one subtract and one division from the critical path. Is that worth it? Only a performance test knows for sure… Feb 15, 2013 at 4:08
• He describes how he got this on page 137-138 with section on "closest point on triangle to point" Aug 9, 2013 at 14:11
• Minor note: there is no argument p to this function.
– Bart
Aug 22, 2014 at 8:05
• Minor implementation note: If all 3 points are on top of each other, you'll get a "divide by 0" error, so be sure to check for that case in the actual code. Feb 8, 2019 at 20:02

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

• Doesn't your proposed solution make assumptions about the third (.z) dimension (specifically, that it doesn't exist)? Sep 29, 2014 at 20:15
• This is the best method here if one's working in 2D. Just a minor improvement: one should compute the reciprocal of the denominator in order to use two multiplications and one division instead of two divisions. May 25, 2016 at 15:20

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

• Could you please write the optimized version of this code (computing barycentric coordinates) for 3D? For a tetrahedron a,b,c,d? Thank you very much Jan 26, 2023 at 10:42

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;


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.

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.

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.

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

}