Finding the centroid of a polygon?

Question

To get the center, I have tried, for each vertex, to add to the total, divide by the number of vertices.

I've also tried to find the topmost, bottommost -> get midpoint... find leftmost, rightmost, find the midpoint.

Both of these did not return the perfect center because I'm relying on the center to scale a polygon.

I want to scale my polygons, so I may put a border around them.

What is the best way to find the centroid of a polygon given that the polygon may be concave, convex and have many many sides of various lengths?

perturbing vertices along the direction from the vertex to the centroid allows you to produce a "rough" but cheap scaling effect. It will work for making an outline or something. A quite simple way to scale in or out by a constant "thickness" is to perturb the vertex by the vector (+/-)A+B where A is unit(next vertex - vertex) and B is unit(prev vertex - vertex). This has the benefit of working on concave polygons. — Steven Lu, Feb 16 '12 at 19:02
Putting a border on a polygon by scaling the polygon is not going to work very well. Imagine one side that goes straight outward from the centroid - no amount of scaling would create a border for it. — Mark Ransom, Dec 07 '12 at 23:55

Emile Cormier · Accepted Answer · 2019-12-13T16:04:56.083

The formula is given here for vertices sorted by their occurance along the polygon's perimeter.

For those having difficulty understanding the sigma notation in those formulas, here is some C++ code showing how to do the computation:

#include <iostream>

struct Point2D
{
    double x;
    double y;
};

Point2D compute2DPolygonCentroid(const Point2D* vertices, int vertexCount)
{
    Point2D centroid = {0, 0};
    double signedArea = 0.0;
    double x0 = 0.0; // Current vertex X
    double y0 = 0.0; // Current vertex Y
    double x1 = 0.0; // Next vertex X
    double y1 = 0.0; // Next vertex Y
    double a = 0.0;  // Partial signed area

    // For all vertices except last
    int i=0;
    for (i=0; i<vertexCount-1; ++i)
    {
        x0 = vertices[i].x;
        y0 = vertices[i].y;
        x1 = vertices[i+1].x;
        y1 = vertices[i+1].y;
        a = x0*y1 - x1*y0;
        signedArea += a;
        centroid.x += (x0 + x1)*a;
        centroid.y += (y0 + y1)*a;
    }

    // Do last vertex separately to avoid performing an expensive
    // modulus operation in each iteration.
    x0 = vertices[i].x;
    y0 = vertices[i].y;
    x1 = vertices[0].x;
    y1 = vertices[0].y;
    a = x0*y1 - x1*y0;
    signedArea += a;
    centroid.x += (x0 + x1)*a;
    centroid.y += (y0 + y1)*a;

    signedArea *= 0.5;
    centroid.x /= (6.0*signedArea);
    centroid.y /= (6.0*signedArea);

    return centroid;
}

int main()
{
    Point2D polygon[] = {{0.0,0.0}, {0.0,10.0}, {10.0,10.0}, {10.0,0.0}};
    size_t vertexCount = sizeof(polygon) / sizeof(polygon[0]);
    Point2D centroid = compute2DPolygonCentroid(polygon, vertexCount);
    std::cout << "Centroid is (" << centroid.x << ", " << centroid.y << ")\n";
}

I've only tested this for a square polygon in the upper-right x/y quadrant.

If you don't mind performing two (potentially expensive) extra modulus operations in each iteration, then you can simplify the previous compute2DPolygonCentroid function to the following:

Point2D compute2DPolygonCentroid(const Point2D* vertices, int vertexCount)
{
    Point2D centroid = {0, 0};
    double signedArea = 0.0;
    double x0 = 0.0; // Current vertex X
    double y0 = 0.0; // Current vertex Y
    double x1 = 0.0; // Next vertex X
    double y1 = 0.0; // Next vertex Y
    double a = 0.0;  // Partial signed area

    // For all vertices
    int i=0;
    for (i=0; i<vertexCount; ++i)
    {
        x0 = vertices[i].x;
        y0 = vertices[i].y;
        x1 = vertices[(i+1) % vertexCount].x;
        y1 = vertices[(i+1) % vertexCount].y;
        a = x0*y1 - x1*y0;
        signedArea += a;
        centroid.x += (x0 + x1)*a;
        centroid.y += (y0 + y1)*a;
    }

    signedArea *= 0.5;
    centroid.x /= (6.0*signedArea);
    centroid.y /= (6.0*signedArea);

    return centroid;
}

I realize this question is old, but instead of the expensive modulus operation, a simple `if(i == nVertexCount - 1)` would have sufficed to handle the wraparound. — SirGuy, Jul 07 '16 at 17:07
@GuyGreer, the modulus thing wasn't my idea, but an edit that someone else proposed. Using a simple `if` like you suggest would certainly work. My answer was intended as a starting point to help folks implement their own centroid algorithm. — Emile Cormier, Jul 07 '16 at 17:25
It's probably worth noting that the order of the set matters. According to the Wikipedia article that was cited, "In these formulas, the vertices are assumed to be numbered in order of their occurrence along the polygon's perimeter". In fact, if you run the code above with the same exact items, but just in a different order, then you'll get different results. — Dave Sexton, Jan 02 '17 at 20:38
Quick heads up, this does not compute the centroid of degenerate polygons, for which it will return {NaN,NaN}. Not a big deal unless you're using the result in other computations and also have ugly input. — nmr, Jun 29 '18 at 19:32
To make this work correctly, you MUST sort the verticies in a correct order. Otherwise the centroid will be in the wrong position. — Desolator, Nov 27 '19 at 23:06
What is the point to multiply `signedArea` by `0.5` and immediately after - by `6.0`? Why not multiply once by `3.0` instead? There are also a lot of redundant initializations! — nikolay, Jul 22 '20 at 01:58
@nikolay It was a translation of the formulas from the Wikipedia article. You'll see the `1/2` and `1/6` factors in those formulas. The compiler may optimize these redundancies away. Feel free to do whatever micro-optimizations you want. :-) — Emile Cormier, Jul 22 '20 at 19:27

score 14 · Answer 2 · edited Feb 26 '17 at 09:21

The centroid can be calculated as the weighted sum of the centroids of the triangles it can be partitioned to.

Here is the C source code for such an algorithm:

/*
    Written by Joseph O'Rourke
    orourke@cs.smith.edu
    October 27, 1995

    Computes the centroid (center of gravity) of an arbitrary
    simple polygon via a weighted sum of signed triangle areas,
    weighted by the centroid of each triangle.
    Reads x,y coordinates from stdin.  
    NB: Assumes points are entered in ccw order!  
    E.g., input for square:
        0   0
        10  0
        10  10
        0   10
    This solves Exercise 12, p.47, of my text,
    Computational Geometry in C.  See the book for an explanation
    of why this works. Follow links from
        http://cs.smith.edu/~orourke/

*/
#include <stdio.h>

#define DIM     2               /* Dimension of points */
typedef int     tPointi[DIM];   /* type integer point */
typedef double  tPointd[DIM];   /* type double point */

#define PMAX    1000            /* Max # of pts in polygon */
typedef tPointi tPolygoni[PMAX];/* type integer polygon */

int     Area2( tPointi a, tPointi b, tPointi c );
void    FindCG( int n, tPolygoni P, tPointd CG );
int     ReadPoints( tPolygoni P );
void    Centroid3( tPointi p1, tPointi p2, tPointi p3, tPointi c );
void    PrintPoint( tPointd p );

int main()
{
    int n;
    tPolygoni   P;
    tPointd CG;

    n = ReadPoints( P );
    FindCG( n, P ,CG);
    printf("The cg is ");
    PrintPoint( CG );
}

/* 
    Returns twice the signed area of the triangle determined by a,b,c,
    positive if a,b,c are oriented ccw, and negative if cw.
*/
int Area2( tPointi a, tPointi b, tPointi c )
{
    return
        (b[0] - a[0]) * (c[1] - a[1]) -
        (c[0] - a[0]) * (b[1] - a[1]);
}

/*      
    Returns the cg in CG.  Computes the weighted sum of
    each triangle's area times its centroid.  Twice area
    and three times centroid is used to avoid division
    until the last moment.
*/
void FindCG( int n, tPolygoni P, tPointd CG )
{
    int     i;
    double  A2, Areasum2 = 0;        /* Partial area sum */    
    tPointi Cent3;

    CG[0] = 0;
    CG[1] = 0;
    for (i = 1; i < n-1; i++) {
        Centroid3( P[0], P[i], P[i+1], Cent3 );
        A2 =  Area2( P[0], P[i], P[i+1]);
        CG[0] += A2 * Cent3[0];
        CG[1] += A2 * Cent3[1];
        Areasum2 += A2;
    }
    CG[0] /= 3 * Areasum2;
    CG[1] /= 3 * Areasum2;
    return;
}

/*
    Returns three times the centroid.  The factor of 3 is
    left in to permit division to be avoided until later.
*/
void Centroid3( tPointi p1, tPointi p2, tPointi p3, tPointi c )
{
    c[0] = p1[0] + p2[0] + p3[0];
    c[1] = p1[1] + p2[1] + p3[1];
    return;
}

void PrintPoint( tPointd p )
{
    int i;

    putchar('(');
    for ( i=0; i<DIM; i++) {
        printf("%f",p[i]);
        if (i != DIM - 1) putchar(',');
    }
    putchar(')');
    putchar('\n');
}

/*
    Reads in the coordinates of the vertices of a polygon from stdin,
    puts them into P, and returns n, the number of vertices.
    The input is assumed to be pairs of whitespace-separated coordinates,
    one pair per line.  The number of points is not part of the input.
*/
int ReadPoints( tPolygoni P )
{
    int n = 0;

    printf("Polygon:\n");
    printf("  i   x   y\n");      
    while ( (n < PMAX) && (scanf("%d %d",&P[n][0],&P[n][1]) != EOF) ) {
        printf("%3d%4d%4d\n", n, P[n][0], P[n][1]);
        ++n;
    }
    if (n < PMAX)
        printf("n = %3d vertices read\n",n);
    else
        printf("Error in ReadPoints:\too many points; max is %d\n", PMAX);
    putchar('\n');

    return  n;
}

There's a polygon centroid article on the CGAFaq (comp.graphics.algorithms FAQ) wiki that explains it.

score 8 · Answer 3 · edited Nov 16 '17 at 17:22

8

boost::geometry::centroid(your_polygon, p);

edited Nov 16 '17 at 17:22

nbro

15,395
32
113
196

answered May 11 '13 at 11:32

Arlen

6,641
4
29
61

score 3 · Answer 4 · answered Sep 15 '20 at 11:44

Here is Emile Cormier's algorithm without duplicated code or expensive modulus operations, best of both worlds:

#include <iostream>

using namespace std;

struct Point2D
{
    double x;
    double y;
};

Point2D compute2DPolygonCentroid(const Point2D* vertices, int vertexCount)
{
    Point2D centroid = {0, 0};
    double signedArea = 0.0;
    double x0 = 0.0; // Current vertex X
    double y0 = 0.0; // Current vertex Y
    double x1 = 0.0; // Next vertex X
    double y1 = 0.0; // Next vertex Y
    double a = 0.0;  // Partial signed area

    int lastdex = vertexCount-1;
    const Point2D* prev = &(vertices[lastdex]);
    const Point2D* next;

    // For all vertices in a loop
    for (int i=0; i<vertexCount; ++i)
    {
        next = &(vertices[i]);
        x0 = prev->x;
        y0 = prev->y;
        x1 = next->x;
        y1 = next->y;
        a = x0*y1 - x1*y0;
        signedArea += a;
        centroid.x += (x0 + x1)*a;
        centroid.y += (y0 + y1)*a;
        prev = next;
    }

    signedArea *= 0.5;
    centroid.x /= (6.0*signedArea);
    centroid.y /= (6.0*signedArea);

    return centroid;
}

int main()
{
    Point2D polygon[] = {{0.0,0.0}, {0.0,10.0}, {10.0,10.0}, {10.0,0.0}};
    size_t vertexCount = sizeof(polygon) / sizeof(polygon[0]);
    Point2D centroid = compute2DPolygonCentroid(polygon, vertexCount);
    std::cout << "Centroid is (" << centroid.x << ", " << centroid.y << ")\n";
}

score 1 · Answer 5 · answered May 08 '10 at 00:40

1

Break it into triangles, find the area and centroid of each, then calculate the average of all the partial centroids using the partial areas as weights. With concavity some of the areas could be negative.

answered May 08 '10 at 00:40

Ben Voigt

277,958
43
419
720

1

Actually it is easier if you break it into right trapezoids using a clever trick. I'll let you discover it. – Alexandru May 08 '10 at 00:47
1

Not sure how it could get any easier. Just a triangle fan where you pick one vertex to stay constant and then take every pair of other adjacent vertices in order. The cross-product then gives the area with the correct sign. – Ben Voigt May 08 '10 at 00:53
The complexity is O(#edges). You need to recreate a right trapezoid for each edge. – Alexandru Sep 15 '15 at 06:46

Finding the centroid of a polygon?

5 Answers5

Linked