Monday, October 20, 2014

Barycentric Coordinates

In the realm of three-dimensional CFD simulations, we tend to use unstructured meshes based on triangles and tetrahedra. If you want to write a program to modify these meshes or even interpolate values, barycentric coordinates can make life much easier.

For simplicity, I'll start with the two-dimensional case of a triangle, but this can be extended easily into three dimensions.  For a triangle, the barycentric coordinate system is comprised of three dependent basis vectors. Each basis vector extends from the midpoint of one side of the triangle to the opposite corner as shown below.


Notice that each corner of the triangle is now represented by a coordinate triple. In fact, any point \((x,y)\) can be represented by a coordinate triple \((\xi_1,\xi_2,\xi_3)\). We also notice that any point outside the triangle will have at least one negative coordinate. This makes barycentric coordinates extremely useful when determining whether a point is inside a triangle. In addition, the coordinates of a point must always add up to 1.

Barycentric coordinates of a point are quite simple to calculate. If we take a point inside a triangle and connect that point to each vertex of the triangle, the result is three subtriangles with areas \(A_1, A_2,\) and \(A_3\) as shown below. If the total area of the triangle is \(A\), then the barycentric coordinates of point \(P\) are simply \(\left(\frac{A_1}{A},\frac{A_2}{A},\frac{A_3}{A}\right)\).


The barycentric coordinates are simply a ratio of the area of each subtriangle to the total area of the triangle.  Now you can see why the coordinates must add up to one.  The overall area of the triangle is the sum of the areas of the three subtriangles. "But wait," you may say, "this won't work for points outside the triangle." Actually, it does work. You may recall that area is a vector quantity, which means there is a direction associated with it. This means that a point outside the triangle will give a negative area for one or more triangles. Let's take a look at this in more detail.

We will define some vectors to help us out. I'll use the notation \(\vec{\mathbf{V}}_{ab}\) to represent a vector from point \(a\) to point \(b\) in the figure above, and similarly for other vectors. Now, let us define some vectors normal to this triangle:\begin{aligned}
\vec{\mathbf{n}}=\vec{\mathbf{V}}_{ab} \times \vec{\mathbf{V}}_{ac} \\
\vec{\mathbf{n}}_1=\vec{\mathbf{V}}_{bc} \times \vec{\mathbf{V}}_{bp} \\
\vec{\mathbf{n}}_2=\vec{\mathbf{V}}_{ca} \times \vec{\mathbf{V}}_{cp} \\
\vec{\mathbf{n}}_3=\vec{\mathbf{V}}_{ab} \times \vec{\mathbf{V}}_{ap} \end{aligned}
Now, you may recall that the area of a triangle can be calculated as half of the cross product of the vectors forming two sides of the triangle. We can thus write the areas as follows:
\begin{aligned}
A=\frac{1}{2}|\vec{\mathbf{n}}| \\
A_1=\frac{1}{2}|\vec{\mathbf{n}}_1| \\
A_2=\frac{1}{2}|\vec{\mathbf{n}}_2| \\
A_3=\frac{1}{2}|\vec{\mathbf{n}}_3| \\
\end{aligned}
These always give positive values for area, but we can check whether the area is really negative by comparing the normal vector of the subtriangle with the normal vector of the overall triangle. We do this with a dot product. If the normal vectors are in the same direction, the result will be positive. Likewise, if the normal vectors are in opposite directions, the dot product will be negative. The way we have defined our normal vectors ensures the correct orientation. For example,
\begin{equation}\frac{\vec{\mathbf{n}}\cdot\vec{\mathbf{n}}_1}{|\vec{\mathbf{n}}||\vec{\mathbf{n}}_1|}=%
\begin{cases}
1 & \text{if P is inside the triangle} \\
 -1 & \text{if P is outside the triangle, across edge b-c} \end{cases}
\end{equation}
It follows that the barycentric coordinates can be calculated by multiplying the area ratio by the corresponding dot-product term.
\begin{aligned}
\xi_1=\frac{\vec{\mathbf{n}}\cdot\vec{\mathbf{n}}_1}{|\vec{\mathbf{n}}|^2} \\
\xi_2=\frac{\vec{\mathbf{n}}\cdot\vec{\mathbf{n}}_2}{|\vec{\mathbf{n}}|^2} \\
\xi_3=\frac{\vec{\mathbf{n}}\cdot\vec{\mathbf{n}}_3}{|\vec{\mathbf{n}}|^2}
\end{aligned}
This works with a triangle in two or three dimensions. We can do something very similar for a tetrahedron. If the barycentric coordinates of a triangle are based on the area of subtriangles, you may have guessed that the barycentric coordinates of a tetrahedron are based on the volume of subtetrahedra, and the coordinates will have four components. The tetrahedron below, formed by vertices \(a,b,c,d\), can be subdivided by inserting a point \(p\).


The volume of a tetrahedron can be calculated with a scalar triple product. Let \(V\) be the total volume of the tetrahedron, \(V_a\) be the volume of the subtetrahedron opposite vertex \(a\), and similarly for the others. Then,
\begin{aligned}
V_a=\frac{1}{6}(\vec{\mathbf{V}}_{bp}\cdot(\vec{\mathbf{V}}_{bd}\times \vec{\mathbf{V}}_{bc})) \\
V_b=\frac{1}{6}(\vec{\mathbf{V}}_{ap}\cdot(\vec{\mathbf{V}}_{ac}\times \vec{\mathbf{V}}_{ad})) \\
V_c=\frac{1}{6}(\vec{\mathbf{V}}_{ap}\cdot(\vec{\mathbf{V}}_{ad}\times \vec{\mathbf{V}}_{ab})) \\
V_d=\frac{1}{6}(\vec{\mathbf{V}}_{ap}\cdot(\vec{\mathbf{V}}_{ab}\times \vec{\mathbf{V}}_{ac})) \\
V=\frac{1}{6}(\vec{\mathbf{V}}_{ab}\cdot(\vec{\mathbf{V}}_{ac}\times \vec{\mathbf{V}}_{ad}))
\end{aligned}
The barycentric coordinates can then be calculated as
\begin{aligned}
\xi_1=\frac{V_a}{|V|} \\
\xi_2=\frac{V_b}{|V|} \\
\xi_3=\frac{V_c}{|V|} \\
\xi_4=\frac{V_d}{|V|}
\end{aligned}

Both of these algorithms are easily implemented in a computer program, and I am constantly finding uses. Here are a few examples of how you can use barycentric coordinates in your programs.

Barycentric Interpolation

Especially if you work with CFD or finite element codes, you will run into triangular and tetrahedron-based meshes.  Sometimes, we want to know the value of some property between nodes on our mesh. We can do this very quickly with barycentric interpolation. If the value at node \(i\) on your element is \(Q_i\), then the value at a point inside that element is given by
\[Q=\sum_{i=1}^n Q_i\xi_i\]

Locating the Element That Contains a Point

As I mentioned earlier, if any of the barycentric coordinate components are negative, then the point lies outside of your element. Even cooler is the fact that the point will lie on the side of the element with the negative component. This enables us to create a semi-intelligent search algorithm. We can take an initial guess at which element the point is in. If it's not the right element, we simply check the element across the side with the negative barycentric component. If we continue this, we will eventually find an element that gives us all positive barycentric components. That element contains the point.

There are many other uses for barycentric coordinates. If you have used barycentric coordinates for something interesting, let us know in the comments.

Source Code

I have written a simple c library that you can use as you desire. Feel free to use, modify, or distribute it as you see fit. All I ask is you give me credit.