[Graphics] Barycentric Coordinates (무게중심 좌표계)

Barycentric Coordinates

Figure 1: barycentric coordinates can be seen as the area of sub-triangles CAP (for u), ABP (for v) and BCP (for w) over the area of the triangle ABC which is why the are also called areal coordinates.

Barycentric coordinates are very useful in many field and are particularly important in CG. They have a few functions and are the key to the next ray-triangle intersection algorithm proposed by Möller-Trumbore that will be studied in the next chapter. How barycentric coordinates can be used in CG will be discussed at the end of this chapter.

Barycentric coordinates can be used to express the position of any point located on the triangle with three scalars. The location of this point includes any position inside the triangle, any position on any of the three edges of the triangles, or any one of the three triangle's vertices themselves. To compute the position of this point using barycentric coordinates we use the following equation (1):

P=uA+vB+wC

where A B and C are the vertices of a triangle and u, v, and w (the barycentric coordinates), three real numbers (scalars) such that u+v+w=1 (barycentric coordinates are normalized). Note that from two of the coordinates we can find the third one: w=1−u−v from which we can establish that u+v≤1 (we will use this simple but important property later on). Equation 1 defines the position of a point P on the plane of the triangle formed by the vertices A, B and C. The point is within the triangle (A, B, C) if 0≤u,v,w≤1. If any one of the coordinates is less than zero or greater than one, the point is outside the triangle. If any of them is zero, P is on one of the lines joining the vertices of the triangle.

You can also simply use two coordinates (lets say u and v) to express the coordinates of P in a two dimensional coordinate system defined by its origin (A) and the edge AB and AC (pretty much like expressing 2D points in the orthogonal two-dimensional coordinate system defined by an x- and y-axis. The only difference in our case is that AB and AC are not necessarily orthogonal and that the origin of this coordinate system is A). Anyway, just to say that you can define a position inside the triangle with the equation P=A+u∗AB+v∗AC (where u≥0 and v≥0 and u+v≤1). This equation can read as, "starts from A, move a little bit in the direction of AB, then a little bit in the direction of AC and you will find P". Now if you develop this equation you can write:

P=A+u(B−A)+v(C−A)=A+uB−uA+vC−vA=(1−u−v)∗A+u∗B+v∗C

This equation is similar to equation 1 excepted that if w=1−u−v you have the following formula:

P=w∗A+u∗B+v∗C instead of P=u∗A+v∗B+w∗C

These two equations are perfectly similar but it's hard to understand why we usually write(1−u−v)∗A+u∗B+v∗C instead of u∗A+v∗B+(1−u−v)∗C without the above explanation.

Barycentric coordinates are also known as areal coordinates. Although not very commonly used, this term indicates that the coordinates u, v and w are proportional to the area of the three sub-triangles defined by P, the point located on the triangle, and the triangle's vertices (A, B, C). These three sub-triangles are denoted ABP, BCP, CAP (see figure 1).

Which leads us to the formulas used to compute the barycentric coordinates:

u=TriangleCAPAreaTriangleABCAreav=TriangleABPAreaTriangleABCAreaw=TriangleBCPAreaTriangleABCArea

In the rest of the lesson we will assume that u makes us move along the edge AB (if u=1 then P=B) and v along the edge AC (if v=1 then P=C). Which is the reason why we will use the area of the sub-triangle CAP to compute u and the area of the sub-triangle ABP to compute v. This is a convention that most people follow in the CG programming community (when we will learn how to use barycentric coordinates to interpolate vertex data you will have a visual example (figure 3) to better understand this).

Figure 2: to compute the area of a triangle we start from the formula used to compute a parallelogram (its base multiplied by its height H). To compute H, we multiply sin(θ) by the length of the vector AC.

Now computing the area of a triangle is trivial. If you duplicate the triangle and mirror it along its longest edge, you get a parallelogram. To compute the area of a parallelogram, simply compute its base, its side and multiply these two numbers together scaled by sin(θ), where θ is the angle subtended by the vectors AB and AC (figure 2). To create the parallelogram we used two triangles, therefore the area of one triangle is half the area of the parallelogram.

With this in hands, it becomes easy to compute u and v (w is computed out of u and v as explained before). 

Trianglearea=||(B−A)||∗||(C−A)||sin(θ)2

To make things easier, we can take advantage of the fact that the area of the sub-triangles ABP, BCP and CAP are proportional to the lengths of the cross product that we computed in the previous chapter to find P, the intersection point. This is one of the cross product properties: the magnitude of the cross product can be interpreted as the area of the parallelogram. Therefore we don't need to explicitly compute the previous formula (which includes an "expansive" sin( θ). We can simply use instead:

Parallelogramarea= ||(B−A)×(C−A)||Trianglearea=Parallelogramarea2

Note than in mathematical terms, the double bars notation (|| ||) means "length of" (lesson 4). In other words we need to compute the length of the vector resulting of the cross product (C-B)x(P-B). We know how to compute the cross product of two vectors and the length of a vector so we have now all we need to compute the barycentric coordinates of the intersection point:

static inline float length(Vector v) { return sqrtf(v.x*v.x + v.y*v.y + v.z*v.z); } // // Compute the intersection of ray with a triangle using the simplest approach // Input: // A, B, C | point | triangle's vertices in CW order // rayOrig | point | ray's origin // rayDir | vector | ray's direction // t | float | distance from the ray origin to the intersection point // u, v | float | barycentric coordinates of intersection point // Return: // return 1 if the ray intersects the triangle, 0 otherwise // bool intersectTriangle(..., float &u, float &v) { // compute plane's normal vector AB, AC; AB = B - A; AC = C - A; // do not normalize vector N = cross(AB, AC); ... // compute triangle area (needed to compute barycentric coord) float triArea = length(N) / 2; // compute u which is the area of sub-tri APC // divided by area of triangle ABC vector AP = P - A; vector APxAC = cross(AP, AC); u = (length(APxAC) / 2) / triArea; // compute v which is the area of sub-tri APB // divided by area of triangle ABC vector APxAB = cross(AP, AB); v = (length(APxAB) / 2) / triArea; ... }

The plane normal should not be normalized because we use the length of the vector to compute the triangle area.

Using Barycentric Coordinates

Barycentric coordinates are most useful in shading. A triangle is a flat surface and we can associate any additional information or data (points, color, vectors, etc.) to each one of its vertices. This information is usually called vertex data. Let's imaging for example that you want vertex A to be red,  vertex B to be green and vertex C to be blue.

Figure 3: barycentric coordinates can be used to interpolate vertex data at the hit point. In this example for example we compute the color at P using vertex color.

If the intersection point coincides with one of the triangle's vertices, then the color of the object at the intersection point is simply the color associated with that vertex. Simple enough. The question is what happens when the ray intersects the triangle anywhere else (either on an edge or inside the triangle)? If the barycentric coordinates are used to compute the position of a point located on the triangle using the triangle vertices, we can interpolate any other data defined at the triangle's vertices (like for example the color) in the exact same way. In other words, barycentric coordinates are used to interpolate vertex data across the triangle's surface (the technique can be applied to any data type, float, color, etc.). This technique is very useful for shading for example to interpolate the normal at the intersection point. Normals of an object can be defined on per face or per vertex basis (we speak of face normal or vertex normal). If they are defined per vertex, we can use this interpolation technique to simulate a smooth shading across the triangle's surface even though, the triangle is "mathematically" flat (the normal at the hit point is a combination of the vertex normals therefore if the vertex normals are different from each other, the result of this interpolation is not constant across the surface of the triangle).

// vertex position point triVertex[3] = {{-3,-3,5}, {0,3,5}, {3,-3,5}}; // vertex data color triColor[3] = {{1,0,0}, {0,1,0}, {0,0,1}}; if (IntersectTriangle(...)) { // compute pixel color // col = w*col0 + u*col1 + v*col2 where w = 1-u-v color objectColor = (1 - u - v) * triColor[0] + u * triColor[1] + v * triColor[2]; }

Barycentric coordinates are also used to compute texture coordinates (we will study this in the lesson on texturing).

Optimizing The Computation Of Barycentric Coordinates

You will notice that in the version of the algorithm we have been describing so far, we use the cross product of AB and AP, and the cross product of CA and CP to compute u and v. But if you look at the code again you will also notice that we have already computed these cross products for the inside-outside test. They have been used to compute the result of whether P is on the right or left side of edge0 (ABxAP) and edge2 (CAxCP). Naturally the first optimisation consists of reusing the result of these values. Also notice that (equation 2):

v=triangleABPareatriangleABCarea=parellogramABParea/2parellogramABCarea/2=parellogramABPareaparellogramABCarea=||AB×AP||||AB×AC||

There is not need to compute the triangle area since the ratio between the triangle ABP and the triangle ABC is the same as the ratio between the parallelogram ABP (which is twice the area of the triangle ABP) and the paralleogram ABC (which is twice the area of the triangle ABC). Therefore we can also avoid the division by 2. The pseudocode becomes:

bool intersectTriangle(..., float &u, float &v) { // compute plane's normal vector AB, AC; AB = B - A; AC = C - A; // do not normalize vector N = cross(AB, AC); ... // compute triangle area (needed to compute barycentric coord) float parallelogramArea = length(N); // no division by 2 // needed for the inside-out test vector AP = P - A; vector APxAC = cross(AP, AC); u = (length(APxAC)) / parallelogramArea; // no division by 2 // needed for the inside-out test vector APxAB = cross(AP, AB); v = (length(APxAB)) / parallelogramArea; // no division by 2 ... }

Finally we can prove that (equation 3):

v=triangleABPareatriangleABCarea=(AB×AP)⋅N(AB×AC)⋅N

First remember that N=AB×AC. Therefore we can rewrite the above above formula as:

v=triangleABPareatriangleABCarea=(AB×AP)⋅NN⋅N

We now need to prove that this equation is true. Remember from lesson 4 that the cross product of two vectors can be seen as:

A×B=cos(θ)||A||||B||

where the angle θ is the angle subtended by the two vectors A and B and where ||A|| and ||B|| is the length of vector A and B. We also know that when the two vector A and B are colinear (they point in the exact same direction), the subtended angle is 0 therefore cos(0) = 1. The result of N⋅N is therefore (the dot product of a vector with itself) cos(0)||N||||N||=||N||||N||, the square of the vector N's length. Now lets look at the numerotor in equation 3. We know from construction (because the points A, B, C and P are coplanar) that the triangles ABP and ABC are coplanar. Therefore the vectors resulting from the cross product of AB×AP and AB×AC are also colinear. We can rewrite the numerator as:

(AB×AP)⋅N=(AB×AP)⋅(AB×AC)

If we assume that AB×AP=A and AB×AC=N=B, we can also write:

(AB×AP)⋅N=cos(0)||A||||B||=||A||||B||

A and B are also colinear in that case because they are constructed from vectors which are coplanar. Finally we can replace our results for the numerator and the denominator in equations 3 and rewrite:

v=triangleABPareatriangleABCarea=||AB×AP||⋅NN⋅N=||A||||B||||B||||B||

We can eliminate one of the ||B||'s which leaves us with:

v=triangleABPareatriangleABCarea=||A||||B||

If we replace A back with AB×AP and B with AB×AC we get:

||A||||B||=||AB×AP||||AB×AC||

which is indeed the formula we used in the first place to compute v (equation 2). Therefore equation 3 is a valid formula to compute v. We already compute the numerator(AB×AP)⋅N when we do the inside-outside test, therefore we can reuse this result in the computation of the barycentric coordinates (a similar technique can be used for u).

In mathematics the formula (A×B)⋅C is called a scalar triple product (because it is the product of three vectors). It can be see as computing the volume of a parallelepiped which sides are defined by the vectors A B and C. It can also be understood as the determinant of the 3x3 matrix having the three vectors as its rows, a property we will be using in the next chapter.

Here is the final version of our intersection routine which includes the optimized method to compute the barycentric coordinates:

bool intersectTriangle(point A, ...) { // compute plane's normal vector AB, AC; AB = B - A; AC = C - A; // DO NOT normalize vector N = Cross(AB, AC); ... // // Step 2: inside-outside test // // edge 1, result needed to compute v vector AP = P - A; vector ABxAP = cross(AB, AP); float v_num; if ((v_num = dot(N, ABxAP)) < 0) return 0; // P is on the left side // edge 2 vector BP = P - B; vector BC = C - B; vector BCxBP = cross(BC, BP); if (Dot(N, BCxBP) < 0) return 0; // P is on the left side // edge 3, needed to compute u vector CP = P - C; // we have computed AC already so we can avoid computing CA by // inverting the vectors in the cross product: // Cross(CA, CP) = cross(CP, AC); vector CAxCP = cross(CP, AC); float u_num; if ((u_num = dot(N, CAxCP)) < 0) return 0; // P is on the left side; // compute barycentric coordinates float denom = dot(N, N); // ABxAC.N where N = ABxAC u = u_num / denom; v = v_num / denom; return 1; // this ray hits the triangle }

Source Code

In this version of the ray-triangle intersection method we have implemented the technique described in chapter 2 to compute the hit point coordinates and the method described in this chapter to compute barycentric coordinates. When a ray hits the triangle, the pixel's color is set with the barycentric coordinates (to help you visualize them).

This code could be optimised even further. We know that the intersection point is not contained by the triangle if u < 0 or u > 1 or v < 0 or v > 1 or (u + v) > 1. Therefore we could first compute u and v and return false as soon as either one of them is lower than 0 or greater than 1 or if (u + v) is greater than 1. Lines 39 to 42 could therefore be entirely skipped. This optimisation is used in the ray-triangle intersection method we present in the next chapter. The IsectData structure was extended to support barycentric coordinates. When an intersection occurs, we set its u and v values (lines 52-53).

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82

struct IsectData { float t; float u, v; }; class Triangle : public Object { public: Triangle(const Matrix44f &o2w, bool singleSided = true) : Object(o2w), isSingledSided(singleSided) { v0 = Vec3f(-1, -1, 0); v1 = Vec3f( 1, -1, 0); v2 = Vec3f( 0, 1, 0); } bool intersect(const Ray<float> &r, IsectData &isectData) const { Vec3f v0v1 = v1 - v0; Vec3f v0v2 = v2 - v0; Vec3f N = cross(v0v1, v0v2); float nDotRay = dot(N, r.dir); if (nDotRay == 0 || (nDotRay > 0 && isSingledSided)) return false; // ray parallel to triangle float d = dot(N, v0); float t = -(dot(N, r.orig) + d) / nDotRay; if (t < 0) return false; // ray behind triangle // inside-out test Vec3f Phit = r(t); // inside-out test edge0 Vec3f v0p = Phit - v0; float v = dot(N, cross(v0v1, v0p)); if (v < 0) return false; // P outside triangle // inside-out test edge1 Vec3f v1p = Phit - v1; Vec3f v1v2 = v2 - v1; float w = dot(N, cross(v1v2, v1p)); if (w < 0) return false; // P outside triangle // inside-out test edge2 Vec3f v2p = Phit - v2; Vec3f v2v0 = v0 - v2; float u = dot(N, cross(v2v0, v2p)); if (u < 0) return false; // P outside triangle float nlen2 = dot(N, N); isectData.t = t; isectData.u = u / nlen2; isectData.v = v / nlen2; return true; } float d; Vec3f v0, v1, v2; bool isSingledSided; }; template<typename T> void trace(const Ray<T> ray, const RenderContext *rc, Spectrum &s) { float tClosest = ray.tmax; Vec3f objectColor; Object *hitObject = NULL; IsectData isectData; for (size_t i = 0; i < rc->objects.size(); ++i) { IsectData isectDataCurrent; if (rc->objects[i]->intersect(ray, isectDataCurrent)) { if (isectDataCurrent.t < tClosest && isectDataCurrent.t > ray.tmin) { isectData = isectDataCurrent; hitObject = rc->objects[i]; objectColor = rc->objects[i]->color; tClosest = isectDataCurrent.t; } } } // an other choice is: if (hitObject != NULL) s.xyz = (tClosest < ray.tmax) ? Vec3f(1 - isectData.u - isectData.v, isectData.u, isectData.v) : rc->backgroundColor; }

Figure 4 shows the result of the program.

Figure 4: a ray traced triangle. Colors were set at the vertices and interpolated across the surface of the triangle using barycentric coordinates.