# 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. – Peter Taylor Feb 12 '12 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. – ashes999 Feb 13 '12 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." – Sean Middleditch Feb 14 '13 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… – John Calsbeek Feb 15 '13 at 4:08
• He describes how he got this on page 137-138 with section on "closest point on triangle to point" – bobobobo Aug 9 '13 at 14:11
• Minor note: there is no argument p to this function. – Bart Aug 22 '14 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. – frodo2975 Feb 8 '19 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)? – Cornstalks Sep 29 '14 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. – rubik May 25 '16 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;
}


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.

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;


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.