Mapping 3d points to 2d and Polygonal Centroids

I ran into a problem while trying to figure out how to get the centroid for an arbitrary convex polygon in 3d. As it turns out internet mathematicians like to be coy so finding a straight answer isn't easy so hopefully someone who needs this will stumble here and think I did a better job.

The math here should be useful to do coordinate conversions from 3D to 2D and back again for other purposes as well.

TLDR: See code at the bottom.

Understanding the problem

A "centroid" is the middle point of the polygon assuming all points are of equal weight. All I want is a centroid for an arbitrary convex polygon rectangles, pentagons, octogons etc. However my coordinates are not 2d, these polygons while flat are floating out in space. So the output is a 3d coordinate.

Looking up the centroid for a polygon you'll likely be pointed to The centroid article on Wikipedia:

This isn't the worst example of Greek alphabet soup, but it's still intimidating and probably not immediately helpful for the budding 3d programmer who's just trying to make stuff work. The second problem is that this only gives it in 2d coordinates. What do we do for 3d? Well according the the accepted answer on Stack Overflow:

Gee, thanks. Well how do we do that? Looking up that gets you a lot half answers and very few good ones.

Anyway we can break down the problem as:

1) Convert 3d coordinates to 2D planar coordinates
2) Calculate the centroid using the above equation
3) Convert back to 3d

Converting co-planar coordinates in 3D to 2D

So we have a bunch of points and we know they are in the same plane. We want to get some 2D coordinates for them. This is a somewhat common operation though you'll see it done in lots of different ways UV coordinates, projection matrix etc. But the generalized version of this is:

1) Make a 2D coordinate basis
2) Map the points

Make a 2d coordinate basis

A basis in this case is a set of vectors that represent what a "step" is the various directions. For 2D we have 2 vectors, you can call them X and Y but because we're doing conversions between a coordinate system that already has an X and Y this might be confusing, we'll call them U and V which is a notation common to texture coordinates. The only rule here is that the vectors of the basis are orthogonal (a change in one will not produce change in the other).

So how do we find a basis? Let's say we have 3 points, the minimum number of points to make up a planar polygon. First we find the normal of the plane. We can do this with the cross product. Given points A,B, and C we make 2 vectors: AB and AC.

function subtractVector(a, b) {
    return [
        a[0] - b[0],
        a[1] - b[1],
        a[2] - b[2]
    ];
}
function crossVector(a, b) {
    return [
        a[1] * b[2] - a[2] * b[1],
        a[2] * b[0] - a[0] * b[2],
        a[0] * b[1] - a[1] * b[0]
    ];
}
function triangleNormal(pointA, pointB, pointC){
    const vector1 = subtractVector(pointC, pointA);
    const vector2 = subtractVector(pointB, pointA);
    return normalizeVector(crossVector(vector1, vector2));
}

The cross product gets us a vector that is orthogonal to 2 vectors so it doesn't really matter what vectors we use as long as they are in the plane we will get a normal. For a complex polygon (more than 3 points) we can just pick any combination of them. However, the cross product is order sensitive. This expect points to be counter-clockwise in order, if not you might get a vector pointing the opposite direction. To be more generic we should normalize the length too.

function normalizeVector(vec){
  const magnitude = Math.sqrt(vec[0]**2 + vec[1]**2 + vec[2]**2);
  return [vec[0] / magnitude, vec[1] / magnitude, vec[2] / magnitude];
}

This works for 3d but you can add terms in the square root for vectors of 4, 5 etc. You take the magnitude which is the length given by square-rooting the sum of squares of each term (Pythagoras's theorem). The we divide each component by that value. What this does is produce a vector of length 1.

Now we have a normal. And that normal is guaranteed to be orthogonal to the first vector AB. Well we can just do that again to get another vector orthogonal to both!

const n = triangleNormal(points[0], points[1], points[2]);
const u = normalizeVector(subtractVector(points[1], points[0])); //ab
const v = normalizeVector(crossVector(u, n));

Note that it doesn't matter which points you pick to get u as long as they are in the plane. The coordinate system might change based on which things you picked but u v and a chosen origin point will let us convert back so it doesn't matter what the coordinate system looks like during calculation, just the fact that it is 2d is enough. u and v should be normalized as well. Now we have our coordinate system: u and v, orthogonal vectors in the plane each with length 1.

Map the points

This part is easy it's the dot product! The dot product measures the similarity between to vectors (eg. orthogonal vectors are 0, vectors in the same direction are 1). It can also be looked at as "projecting" one vector on to another which is exactly our case. We project each point onto the u and v basis vectors.

const n = triangleNormal(points[0], points[1], points[2]);
const u = normalizeVector(subtractVector(points[1], points[0]));
const v = normalizeVector(crossVector(u, n));
const p0 = points[0];

const mappedPoints = points.map(p => [dotVector(subtractVector(p, p0),u), dotVector(subtractVector(p, p0),v)]);

export function dotVector(a, b) {
    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}

mappedPoints contains our 3d coordinates in the 2d UV space. The subtraction in the dotVector is what sets up the origin point. All points will use p0 as the origin so we're effectively calculating them all relative to that origin point. This is needed because p0 exists in both the plane and 3d space and is the extra bit of information we need to convert back to 3d. Again it doesn't matter what the point is (which is why I took the first) it can be anything as long as it exists in both spaces.

Calculate the centroid

This is mostly just transcribing moon language from Wikipedia. There are 3 steps:

1) Calculate area
2) Calculate CX
3) Calculate CY

Calculate the Area

You can reference the equation from wikipedia which gives the shoe-lace formula. I have no idea how it's derived but thankfully the page provides some examples so we can actually test our implementations! Anyway, aside from being a very clever way to do things on paper it's just a loop:

export function polyArea(points){
    let sum = 0;
    for(let i = 0; i < points.length; i++){
        const nextI = (i + 1) % points.length;
        sum += points[i][0] * points[nextI][1] - points[nextI][0] * points[i][1];
    }
    return Math.abs(sum) / 2;
}

We need to get the i value in one column, the next i value of the other column, add them and switch directions and subtract that total. In this case we do the subtraction in the same step. nextI here ensures that we wrap around as the last i in one column corresponds to the first i in the other. In the end we halve the absolute value. Note that the absolute value is useful if this is a generic function because area is always positive, but it actually isn't necessary for the centroid calculation to be correct.

Calculate the 2d centroid coordinates

Again it's just implementing the equation making sure that the points wrap around:

export function polyCentroid2d(points){
    const area = polyArea(points);

    let sumX = 0;
    let sumY = 0;
    for (let i = 0; i < points.length; i++) {
        const nextI = (i + 1) % points.length;
        const x0 = points[i][0];
        const x1 = points[nextI][0];
        const y0 = points[i][1];
        const y1 = points[nextI][1];

        const doubleArea = (x0 * y1) - (x1 * y0);
        sumX += (x0 + x1) * doubleArea;
        sumY += (y0 + y1) * doubleArea;
    }

    const cx = sumX / (6 * area);
    const cy = sumY / (6 * area);

    return [cx, cy];
}   

The thing that sucks about compact equations is that you often don't know what to call something. Would you have figured out that (x0 * y1) - (x1 * y0) is an area calculation from the original equation? Probably not. But we can recognize this from the polyArea function, the same terms in the same loop. We divided that sum by 2 because it's double the area. Not intuitive at all. But if you don't need polyArea for anything else you can fold that logic into the loop here too which is how it's done in the final code.

Converting 2D planar coordinates back to 3D

To go back we need some sort of inverse to the mapping operation.

const [cu, cv] = polyCentroid2d(mappedPoints);

const x = dotVector([p0[0], u[0], v[0]], [1, cu, cv]),
const y = dotVector([p0[1], u[1], v[1]], [1, cu, cv]),
const z = dotVector([p0[2], u[2], v[2]], [1, cu, cv])

This time we use our basis vectors to map backwards. However since we lost a dimension worth of info in the conversion we need to use the origin point mentioned earlier to recover that and get the final X, Y and Z coordinate for the centroid.

The code

I combined the area calculation with the centroid calculation to make it a little more compact: