# Find the Centroid of a Polygon

From Wikipedia:

The centroid of a non-self-intersecting closed polygon defined by n vertices (x0,y0), (x1,y1), ..., (xn﻿−﻿1,yn−1) is the point (Cx, Cy), where

and where A is the polygon's signed area,

In these formulas, the vertices are assumed to be numbered in order of their occurrence along the polygon's perimeter. Furthermore, the vertex ( xn, yn ) is assumed to be the same as ( x0, y0 ), meaning i + 1 on the last case must loop around to i = 0. Note that if the points are numbered in clockwise order the area A, computed as above, will have a negative sign; but the centroid coordinates will be correct even in this case.

• Given a list of vertices in order (either clockwise, or counter-clockwise), find the centroid of the non-self-intersecting closed polygon represented by the vertices.
• If it helps, you may assume input to be only CW, or only CCW. Say so in your answer if you require this.
• The coordinates are not required to be integers, and may contain negative numbers.
• Input will always be valid and contain at least three vertices.
• Inputs only need to be handled that fit in your language's native floating point data type.
• You may assume that input numbers will always contain a decimal point.
• You may assume that input integers end in . or .0.
• You may use complex numbers for input.
• Output should be accurate to the nearest thousandth.

### Examples

[(0.,0.), (1.,0.), (1.,1.), (0.,1.)]        -> (0.5, 0.5)
[(-15.21,0.8), (10.1,-0.3), (-0.07,23.55)]  -> -1.727 8.017
[(-39.00,-55.94), (-56.08,-4.73), (-72.64,12.12), (-31.04,53.58), (-30.36,28.29), (17.96,59.17), (0.00,0.00), (10.00,0.00), (20.00,0.00), (148.63,114.32), (8.06,-41.04), (-41.25,34.43)]   -> 5.80104769975, 15.0673812762


Too see each polygon on a coordinate plane, paste the coordinates without the square brackets in the "Edit" menu of this page.

I confirmed my results using this Polygon Centroid Point Calculator, which is awful. I couldn't find one that you can input all vertices at once, or that didn't try to erase your - sign when you type it first. I'll post my Python solution for your use after people have had a chance to answer.

• The much simpler technique of averaging all x's and y's works for the first two sets, but not the third. I wonder what makes the difference... Oct 18 '16 at 22:23
• @ETHproductions The third polygon is not convex. Oct 19 '16 at 0:35
• @ETHproductions If you approximate a circle with a polygon, you can move the average point arbitrary close to a point on the circle by using more points close to that point, while nearly not effecting the centroid and keeping the polygon convex. Oct 19 '16 at 1:04
• @ETHproductions Actually convexity is not the reason. Averaging all xs and ys puts all the weight in the vertices instead of distributed over the body. The first one happen to work because it is regular, so both methods end up at the symmetry center. The second one works because for triangles both methods lead to the same point. Oct 19 '16 at 6:14
• Can we use complex numbers for I/O?
– xnor
Oct 19 '16 at 7:03

# Mathematica, 23 bytes

RegionCentroid@*Polygon


Take THAT, Jelly!

Edit: One does not simply beat Jelly...

# Explanation

Polygon


Generate a polygon with vertices at the points specified.

RegionCentroid


Find the centroid of the polygon.

• Well you beat me, but there's probably a shorter way than what I have, I don't have a complete understanding of Jelly yet Oct 19 '16 at 0:24
• @miles aw... :( Oct 19 '16 at 0:25

# Jelly, 25242221 18 bytes

S×3÷@×"
ṙ-żµÆḊçS€S


Applies the formula shown in the problem.

Saved 3 bytes with help from @Jonathan Allan.

## Explanation

S×3÷@×"  Helper link. Input: determinants on LHS, sum of pairs on RHS
S        Sum the determinants
×3      Multiply by 3
×"  Vectorized multiply between determinants and sums
÷@    Divide that by the determinant sum multipled by 3 and return

ṙ-żµÆḊçS€S  Main link. Input: 2d list of points
ṙ-          Rotate the list of points by 1 to the right
ż         Interleave those with the original points
This creates all overlapping slices of length 2
ÆḊ      Get the determinant of each slice
S€   Get the sum of each slice (sum of pairs of points)
S  Sum and return

• You can replace ṁL‘$ṡ2 with ṙ1ż@ or żṙ1$ Oct 19 '16 at 1:22
• @JonathanAllan Thanks, also I can rotate by ṙ-ż to avoid the swap and save another byte Oct 19 '16 at 1:29
• Oh yes of course! Oct 19 '16 at 1:31

# J, 29 bytes

2+/@(+/$$*%3*1#.])-/ .*$$],{.


Applies the formula shown in the problem.

## Usage

   f =: 2+/@(+/$$*%3*1#.])-/ .*$$],{.
f 0 0 , 1 0 , 1 1 ,: 0 1
0.5 0.5
f _15.21 0.8 , 10.1 _0.3 ,: _0.07 23.55
_1.72667 8.01667
f _39 _55.94 , _56.08 _4.73 , _72.64 12.12 , _31.04 53.58 , _30.36 28.29 , 17.96 59.17 , 0 0 , 10 0 , 20 0 , 148.63 114.32 , 8.06 _41.04 ,: _41.25 34.43
5.80105 15.0674


## Explanation

2+/@(+/$$*%3*1#.])-/ .*$$],{.  Input: 2d array of points P [[x1 y1] [x2 y2] ...]
]     Get P
,    Join, makes the end cycle back to the front
2                              The constant 2
2                      \       For each pair of points
-/ .*        Take the determinant
2    +/\                       Sum each pair of points
*                     Multiply the sum of each pair by its determinant
%                    Divide each by
1#.]              The sum of the determinants
3*                  Multiplied by 3
+/@                           Sum and return


## Maxima, 124 118 116 112 106 byte

f(l):=(l:endcons(l[1],l),l:sum([3,l[i-1]+l[i]]*determinant(matrix(l[i-1],l[i])),i,2,length(l)),l[2]/l[1]);


I'm not experienced with Maxima, so any hints are welcome.

Usage:

(%i6) f([[-15.21,0.8], [10.1,-0.3], [-0.07,23.55]]);
(%o6)              [- 1.726666666666668, 8.016666666666668]


## Racket 420 bytes

(let*((lr list-ref)(getx(lambda(i)(lr(lr l i)0)))(gety(lambda(i)(lr(lr l i)1)))(n(length l))(j(λ(i)(if(= i(sub1 n))0(add1 i))))
(A(/(for/sum((i n))(-(*(getx i)(gety(j i)))(*(getx(j i))(gety i))))2))
(cx(/(for/sum((i n))(*(+(getx i)(getx(j i)))(-(*(getx i)(gety(j i)))(*(getx(j i))(gety i)))))(* 6 A)))
(cy(/(for/sum((i n))(*(+(gety i)(gety(j i)))(-(*(getx i)(gety(j i)))(*(getx(j i))(gety i)))))(* 6 A))))
(list cx cy))


Ungolfed:

(define(f l)
(let* ((lr list-ref)
(getx (lambda(i)(lr (lr l i)0)))
(gety (lambda(i)(lr (lr l i)1)))
(n (length l))
(j (lambda(i) (if (= i (sub1 n)) 0 (add1 i))))
(A (/(for/sum ((i n))
(-(* (getx i) (gety (j i)))
(* (getx (j i)) (gety i))))
2))
(cx (/(for/sum ((i n))
(*(+(getx i)(getx (j i)))
(-(*(getx i)(gety (j i)))
(*(getx (j i))(gety i)))))
(* 6 A)))
(cy (/(for/sum ((i n))
(*(+(gety i)(gety (j i)))
(-(*(getx i)(gety (j i)))
(*(getx (j i))(gety i)))))
(* 6 A))))
(list cx cy)))


Testing:

(f '[(-15.21 0.8)  (10.1 -0.3)  (-0.07 23.55)] )
(f '[(-39.00 -55.94)  (-56.08 -4.73)  (-72.64 12.12)  (-31.04 53.58)
(-30.36 28.29)  (17.96 59.17)  (0.00 0.00)  (10.00 0.00)
(20.00 0.00) (148.63 114.32)  (8.06 -41.04)  (-41.25 34.43)])


Output:

'(-1.7266666666666677 8.01666666666667)
'(5.8010476997538465 15.067381276150996)


# R, 129 127 bytes

function(l){s=sapply;x=s(l,[,1);y=s(l,[,2);X=c(x[-1],x[1]);Y=c(y[-1],y[1]);p=x*Y-X*y;c(sum((x+X)*p),sum((y+Y)*p))/sum(p)/3}


Unnamed function that take an R-list of tuples as input. The named equivalent can be called using e.g.:

f(list(c(-15.21,0.8),c(10.1,-0.3),c(-0.07,23.55)))


Ungolfed and explained

f=function(l){s=sapply;                           # Alias for sapply
x=s(l,[,1);                       # Split list of tuples into vector of first elements
y=s(l,[,2);                       # =||= but for second element
X=c(x[-1],x[1]);                    # Generate a vector for x(i+1)
Y=c(y[-1],y[1]);                    # Generate a vector for y(i+1)
p=x*Y-X*y;                          # Calculate the outer product used in both A, Cx and Cy
c(sum((x+X)*p),sum((y+Y)*p))/sum(p)/3    # See post for explanation
}


The final step (c(sum((x+X)*p),sum((y+Y)*p))/sum(p)*2/6) is a vectorized way of calculating both Cx and Cy. The sum in the formulas for Cx and Cy are stored in a vector and consequently divided by the "sum in A" *2/6. E.g.:

(SUMinCx, SUMinCy) / SUMinA / 3


, and then implicitly printed.

Try it on R-fiddle

• *2/6 could probably be /3? Oct 19 '16 at 19:54
• @mbomb007 It's so painstakingly obvious, I guess I got caught up in golfing the other part. /shrug Oct 20 '16 at 6:55
• Elegant, I like your use of sapply to deal with those lists! There could be scope for golfing here, I'm not sure how flexible the allowable input is. If you're allowed to input just a sequence of coordinates, like c(-15.21,0.8,10.1,-0.3,-0.07,23.55), then you can save 17 bytes by replacing the first lines of your function with y=l[s<-seq(2,sum(1|l),2)];x=l[-s];. That is, setting y to be every even-indexed element of l, and x to be every odd-indexed element. Oct 20 '16 at 7:13
• Even better, though, would be if we can input a matrix (or array), like matrix(c(-15.21,0.8,10.1,-0.3,-0.07,23.55),2), as then the beginning of your function can be x=l[1,];y=l[2,];, which saves 35 bytes. (The input matrix could be transposed, in which case x=l[,1];y=l[,2];.) Of course, the easiest solution of all is if the x and y points are just input as separate vectors, function(x,y), but I don't think that's allowed... Oct 20 '16 at 7:13
• @rturnbull I asked OP in comments and he specifically wanted a list of tuples (very inconvenient in R of course) so I don't think the matrix approach is allowed. And even if it was, the input would have to be the vector part (i.e. c(...)) and the matrix conversion would have to be done inside the function. Oct 20 '16 at 7:33

# Python, 156 127 bytes

def f(p):n=len(p);p=p+p[:1];i=s=0;exec'd=(p[i].conjugate()*p[i+1]).imag;s+=d;p[i]=(p[i]+p[i+1])*d;i+=1;'*n;print sum(p[:n])/s/3


Ungolfed:

def f(points):
n = len(points)
points = points + [points[0]]
determinantSum = 0
for i in range(n):
determinant = (points[i].conjugate() * points[i+1]).imag
determinantSum += determinant
points[i] = (points[i] + points[i+1]) * determinant
print sum(points[:n]) / determinantSum / 3


Ideone it.

This takes each pair of points [x, y] as a complex number x + y*j, and outputs the resulting centroid as a complex number in the same format.

For the pair of points [a, b] and [c, d], the value a*d - b*c which is needed for each pair of points can be computed from the determinant of the matrix

| a b |
| c d |


Using complex arithmetic, the complex values a + b*j and c + d*j can be used as

conjugate(a + b*j) * (c + d*j)
(a - b*j) * (c + d*j)
(a*c + b*d) + (a*d - b*c)*j


Notice that the imaginary part is equivalent to the determinant. Also, using complex values allows the points to be easily summed component-wise in the other operations.

# R + sp (46 bytes)

Assumes sp package is installed (https://cran.r-project.org/web/packages/sp/)

Takes a list of vertices, (for example list(c(0.,0.), c(1.,0.), c(1.,1.), c(0.,1.)))

Takes advantage of the fact that the "labpt" of a Polygon is the centroid.

function(l)sp::Polygon(do.call(rbind,l))@labpt


# JavaScript (ES6), 102

Straight implementation of the formula

l=>[...l,l[0]].map(([x,y],i)=>(i?(a+=w=t*y-x*u,X+=(t+x)*w,Y+=(u+y)*w):X=Y=a=0,t=x,u=y))&&[X/3/a,Y/3/a]


Test

f=
l=>[...l,l[0]].map(([x,y],i)=>(i?(a+=w=t*y-x*u,X+=(t+x)*w,Y+=(u+y)*w):X=Y=a=0,t=x,u=y))&&[X/3/a,Y/3/a]

function go()
{
var c=[],cx,cy;
// build coordinates array
I.value.match(/-?[\d.]+/g).map((v,i)=>i&1?t[1]=+v:c.push(t=[+v]));
console.log(c+''),
[cx,cy]=f(c);
O.textContent='CX:'+cx+' CY:'+cy;
// try to display the polygon
var mx=Math.max(...c.map(v=>v[0])),
nx=Math.min(...c.map(v=>v[0])),
my=Math.max(...c.map(v=>v[1])),
ny=Math.min(...c.map(v=>v[1])),
dx=mx-nx, dy=my-ny,
ctx=C.getContext("2d"),
cw=C.width, ch=C.height,
fx=(mx-nx)/cw, fy=(my-ny)/ch, fs=Math.max(fx,fy)
C.width=cw
ctx.setTransform(1,0,0,1,0,0);
ctx.beginPath();
c.forEach(([x,y],i)=>ctx.lineTo((x-nx)/fs,(y-ny)/fs));
ctx.closePath();
ctx.stroke();
ctx.fillStyle='#ff0000';
ctx.fillRect((cx-nx)/fs-2,(cy-ny)/fs-2,5,5);
}
go()
#I { width:90% }
#C { width:90%; height:200px;}
<input id=I value='[[-15.21,0.8], [10.1,-0.3], [-0.07,23.55]]'>
<button onclick='go()'>GO</button>
<pre id=O></pre>
<canvas id=C></canvas>

# Python 2, 153 bytes

Uses no complex numbers.

P=input()
A=x=y=0;n=len(P)
for i in range(n):m=-~i%n;a=P[i][0];b=P[i][1];c=P[m][0];d=P[m][1];t=a*d-b*c;A+=t;x+=t*(a+c);y+=t*(b+d)
k=1/(3*A);print x*k,y*k


Try it online

### Ungolfed:

def centroid(P):
A=x=y=0
n=len(P)
for i in range(n):
m=-~i%n
x0=P[i][0];y0=P[i][1]
x1=P[m][0];y1=P[m][1]
t = x0*y1 - y0*x1
A += t/2.
x += t * (x0 + x1)
y += t * (y0 + y1)
k = 1/(6*A)
x *= k
y *= k
return x,y


# Actually, 4540 39 bytes

This uses an algorithm similar to miles' Jelly answer. There is a shorter way to calculate determinants using a dot product, but there's currently a bug with Actually's dot product where it won't work with lists of floats. Golfing suggestions welcome. Try it online!

;\Z♂#;i¥M@i│N@F*)F@N*-M;Σ3*)♀*┬♂Σ♀/


Ungolfing

         Implicit input pts.
;\       Duplicate pts, rotate right.
Z        Zip rot_pts and pts together.
♂#       Convert the iterables inside the zip to lists
(currently necessary due to a bug with duplicate)
;        Duplicate the zip.
...M   Get the sum each pair of points in the zip.
i        Flatten the pair to the stack.
¥        Pairwise add the two coordinate vectors.
@        Swap with the other zip.
...M   Get the determinants of the zip.
i│       Flatten to stack and duplicate entire stack.
Stack: [a,b], [c,d], [a,b], [c,d]
N@F*)    Push b*c and move it to BOS.
F@N*     Push a*d.
-        Get a*d-b*c.
;Σ3*)    Push 3 * sum(determinants) and move it to BOS.
♀*       Vector multiply the determinants and the sums.
┬        Transpose the coordinate pairs in the vector.
♂Σ       Sum the x's, then the y's.
♀/       Divide the x and y of this last coordinate pair by 3*sum(determinants).
Implicit return.


A shorter, non-competitive version

This is another 24-byte version that uses complex numbers. It's non-competitive because it relies on bug fixes that postdate this challenge. Try it online!

;\│¥)Ziá*╫@XM;Σ3*)♀*Σ/


Ungolfing

         Implicit input a list of complex numbers, pts.
;\       Duplicate pts, rotate right.
│        Duplicate stack. Stack: rot_pts, pts, rot_pts, pts.
¥)       Pairwise sum the two lists of points together and rotate to BOS.
Z        Zip rot_pts and pts together.
...M   Map the following function over the zipped points to get our determinants.
i        Flatten the list of [a+b*i, c+d*i].
á        Push the complex conjugate of a+bi, i.e. a-b*i.
*        Multiply a-b*i by c+d*i, getting (a*c+b*d)+(a*d-b*c)*i.
Our determinant is the imaginary part of this result.
╫@X      Push Re(z), Im(z) to the stack, and immediately discard Re(z).
This map returns a list of these determinants.
;        Duplicate list_determinants.
Σ3*)     Push 3 * sum(list_determinants) and rotate that to BOS.
♀*Σ      Pairwise multiply the sums of pairs of points and the determinants and sum.
/        Divide that sum by 3*sum(list_determinants).
Implicit return.


# C++14, 241 bytes

struct P{float x;float y;};
#define S(N,T)auto N(P){return 0;}auto N(P a,P b,auto...V){return(T)*(a.x*b.y-b.x*a.y)+N(b,V...);}
S(A,1)S(X,a.x+b.x)S(Y,a.y+b.y)auto f(auto q,auto...p){auto a=A(q,p...,q)*3;return P{X(q,p...,q)/a,Y(q,p...,q)/a};}


Output is the helper struct P,

Ungolfed:

 //helper struct
struct P{float x;float y;};

//Area, Cx and Cy are quite similar
#define S(N,T)\  //N is the function name, T is the term in the sum
auto N(P){return 0;} \   //end of recursion for only 1 element
auto N(P a,P b,auto...V){ \ //extract the first two elements
return (T)*(a.x*b.y-b.x*a.y) //compute with a and b
+ N(b,V...); \        //recursion without first element
}

//instantiate the 3 formulas
S(A,1)
S(X,a.x+b.x)
S(Y,a.y+b.y)

auto f(auto q,auto...p){
auto a=A(q,p...,q)*3; //q,p...,q appends the first element to the end
return P{X(q,p...,q)/a,Y(q,p...,q)/a};
}


Usage:

f(P{0.,0.}, P{1.,0.}, P{1.,1.}, P{0.,1.})
f(P{-15.21,0.8}, P{10.1,-0.3}, P{-0.07,23.55})


## Clojure, 177156 143 bytes

Update: Instead of a callback I'm using [a b c d 1] as a function and the argument is just a list of indexes to this vector. 1 is used as a sentinel value when calculating A.

Update 2: Not precalculating A at let, using (rest(cycle %)) to get input vectors offset by one.

#(let[F(fn[I](apply +(map(fn[[a b][c d]](*(apply +(map[a b c d 1]I))(-(* a d)(* c b))))%(rest(cycle %)))))](for[i[[0 2][1 3]]](/(F i)(F[4])3)))


Original version:

#(let[F(fn[L](apply +(map(fn[[a b][c d]](*(L[a b c d])(-(* a d)(* c b))))%(conj(subvec % 1)(% 0)))))A(*(F(fn[& l]1))3)](map F[(fn[v](/(+(v 0)(v 2))A))(fn[v](/(+(v 1)(v 3))A))]))


At less golfed stage:

(def f (fn[v](let[F (fn[l](apply +(map
(fn[[a b][c d]](*(l a b c d)(-(* a d)(* c b))))
v
(conj(subvec v 1)(v 0)))))
A (* (F(fn[& l] 1)) 3)]
[(F (fn[a b c d](/(+ a c)A)))
(F (fn[a b c d](/(+ b d)A)))])))


Creates a helper function F which implements the summation with any callback l. For A the callback returns constantly 1 whereas X and Y coordinates have their own functions. (conj(subvec v 1)(v 0)) drops the first element and appends to the end, this way it is easy to keep track of x_i and x_(i+1). Maybe there is still some repetition to be eliminated, especially at the last (map F[....