centroid of a triangle with extension to polygonsby bob jantzen inspired by klaus volpert and the city manager of the Philadelphia
December 2, 2004As long as students can write the equations of straight lines between pairs of points, and can set up the double integrals of the constant 1 and the two coordinates LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJw== and LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJw== over a triangular region in the plane, they can derive the formula for the area and centroid of this triangular region using a computer algebra system, which they could probably never do by hand. In principle nothing prevents this since only integration of linear and quadratic polynomials is involved, but it is a bit like multiplying 3 digit numbers in your head. A real misuse of the brain. And the factoring of the result which is necessary for simplifying to the very simple final formula is not something that can easily be done by hand (I am guessing), even by an experienced mathematician, without resorting to tricks.
The remarkably simple result is that the coordinates of the centroid are just the averages of the coordinates of the three vertices of the triangle, and it would be nice if one could provide a simple geometric proof of this fact.This result then enables one to easily compute the centroid of any polygonal region of the plane by triangulation and a superposition principle: once a point is chosen for the triangulation, one can average the centroids of the triangles which make up the polygon weighted by their fractional areas.restart; with(plots): with(linalg):centroid of a triangular region in the planeSuppose we are given the three vertices of a triangle in the plane. We then write the equations for the line segments connecting themY12:=y1+(y2-y1)/(x2-x1)*(x-x1);
Y13:=y1+(y3-y1)/(x3-x1)*(x-x1);
Y23:=y3+(y3-y2)/(x3-x2)*(x-x3);centroid diagramwith(plots):plot([[1,1],[2,3],[4,2],[1,1]],x=0..4.5,y=0...3.5,color=red,thickness=2): triplot:=%:
display(triplot,
plot([[1,1],[1,0]],color=grey),
plot([[2,3],[2,0]],color=grey),
plot([[4,2],[4,0]],color=grey),
plot([[1/3*(1+2+4),1/3*(1+2+3)]],style=point,color=blue),
textplot([2.7,2,`(X,Y)`]),
textplot([0.6,1,`(x1,y1)`]),
textplot([2,3.2,`(x2,y2)`]),
textplot([4.4,2,`(x3,y3)`])
); trianglediagram:=%:With no loss of generality, we assume the configuration shown in the diagram to set up the iterated integrals for the area A and center of gravity coordinates (X,Y) of the triangle:display(trianglediagram);We integrate 1 over the triangle to get the area A:Int(Int(1,y=Y13..Y12),x=x1..x2)
+Int(Int(1,y=Y13..Y23),x=x2..x3);
value(%):
expand(%):
A:=factor(%);We integrate LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJw== over the triangle:Int(Int(x,y=Y13..Y12),x=x1..x2)
+Int(Int(x,y=Y13..Y23),x=x2..x3);
value(%):
expand(%):
factor(%):
AX:=%;Then we integrate LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2JVEieUYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJw== over the triangle:Int(Int(y,y=Y13..Y12),x=x1..x2)
+Int(Int(y,y=Y13..Y23),x=x2..x3);
value(%):
expand(%):
factor(%):
AY:=%;And divide both by the area to get the center of gravity of the triangle, whose coordinates are just the averages of the endpoint coordinates:X:=simplify(AX/A);
Y:=simplify(AY/A);
principle of superposition and oriented areasThe coordinates of the centroid are quotients whose numerator is the area times the centroid coordinate, and the denominator is the area. If the region of integration in the numerator is broken up into a sum of integrals over a family of subdivisions of the whole, then each contribution separately represents the area of the subregion times the centroid coordinate of the subregion. Thus the quotient represents the weighted average of the centroid coordinates of the subregions where the weight is just the fractional area of the total:(AX[1]+AX[2])/(A[1]+A[2])
=A[1]/(A[1]+A[2])*X[1]+A[2]/(A[1]+A[2])*X[2];Thus if we take any convex polygon and triangulate it, the centroid can be found by this weighted average of the centroids of the individual triangles. Not only that but one can easily see that for any polygon, by using oriented areas, one can achieve the same result.The area of a triangle is equal to half the third component crossproduct of two sides of the triangle (as difference vectors, tails at the common vertex) emanating from a given vertex provided that the first difference vector rotates counterclockwise into the second difference vector so that the right hand rule gives a positive sign for the third component of the crossproduct. If we want to list the sides in clockwise order, we would have to reverse the sign of the crossproduct.A;with(linalg):crossprod([x1,y1,0],[x2,y2,0]);
crossprod([x3-x1,y3-y1,0],[x2-x1,y2-y1,0])[3];
simplify(%/2-A);
centroid of a convex polygon in the plane by triangulationSo we now list the points in clockwise order around a polygon and then choose a point in the interior for the center of the spokes out to these points which divide the polygon into a family of triangles.polygon diagramWe include the first point at the end so we can easily plot them and iterate the triangle pairings:points:=[[5,6],[7,3],[7,1],[4,0],[2,1],[1,2],[2,5],[5,6]];
nops(%)-1;"CenterPointT" for the spokes:cpt:=[4,3];pointsspokes:=[[5,6],cpt,[7,3],cpt,[7,1],cpt,[4,0],cpt,[2,1],cpt,[1,2],cpt,[2,5],cpt,[5,6]]:centerpt:=plot([cpt],style=point,color=blue):
perimeterpoints:=plot(points, color=black, style=point):
perimetersides:=plot(points,color=red,thickness=2):
spokes:=plot(pointsspokes,color=grey):
display(perimeterpoints,perimetersides,centerpt,spokes);
polygondiagram:=%:display(polygondiagram);To use the cross-product we must add a third zero component to all the 2 entry difference vectors of the points in the plane:imbed3:=proc(P)
[P[1],P[2],0];
end proc;The following oriented triangle area function gives the area of the triangle if the three points are in clockwise order around the perimeter of the triangle, the negative of the area if they are in counterclockwise order.OArea:=proc(P,Q,C)
-1/2*crossprod(imbed3(P)-imbed3(C),imbed3(Q)-imbed3(C))[3];
end proc;We test it and calculate the total area:OArea(points[1],points[2],cpt);
OArea(points[2],points[1],cpt);
PArea:=add(OArea(points[i],points[i+1],cpt),i=1..7);We then introduce the triangle centroid and oriented-area weighted centroid functions.Centroid:=proc(P,Q,R)
[1/3*(P[1]+Q[1]+R[1]),1/3*(P[2]+Q[2]+R[2])];
end proc;ACentroid:=proc(P,Q,R)
[1/3*(P[1]+Q[1]+R[1])*OArea(P,Q,R),1/3*(P[2]+Q[2]+R[2])*OArea(P,Q,R)];
end proc;We test it:Centroid(points[1],points[2],cpt);
ACentroid(points[1],points[2],cpt);
OArea(points[1],points[2],cpt);Now we compute the centroid of the polygon, adding up all the oriented-area weighted triangle contributions and dividing by the total area:PCentroid:=add(ACentroid(points[i],points[i+1],cpt),i=1..7)
/PArea;
Next we add the individual centroids to the diagram together with the centroid of the polygon, which is a black cross.plot([PCentroid],style=point,color=black,symbol=CROSS):
PCentroidplot:=%:for i from 1 to 7 do
C[i]:=Centroid(points[i],points[i+1],cpt):
end do:
Centroids:=[seq(C[i],i=1..7)];
plot(Centroids,style=point,color=blue): Centroidsplot:=%:display(polygondiagram,Centroidsplot,PCentroidplot);By chance we picked our spoke center point very close to the actual centroid!extension to general polygons?Now we check that we get the same result even if our center point for the spokes does not even lie inside the polygon:cpt:=[0,0]:
pointsspokes:=[[5,6],cpt,[7,3],cpt,[7,1],cpt,[4,0],cpt,[2,1],cpt,[1,2],cpt,[2,5],cpt,[5,6]]:
for i from 1 to 7 do
C[i]:=Centroid(points[i],points[i+1],cpt):
end do:
Centroids:=[seq(C[i],i=1..7)];
plot(Centroids,style=point,color=blue): Centroidsplot:=%:
centerpt:=plot([cpt],style=point,color=blue):
perimeterpoints:=plot(points, color=black, style=point):
perimetersides:=plot(points,color=red,thickness=2):
spokes:=plot(pointsspokes,color=grey):
display(perimeterpoints,perimetersides,centerpt,spokes):
polygondiagram:=%:
display(polygondiagram,Centroidsplot,PCentroidplot);Now we get both the correct area and centroid for the whole:add(OArea(points[i],points[i+1],cpt),i=1..7);
add(ACentroid(points[i],points[i+1],cpt),i=1..7)/PArea;
One gets a cancellation of the contributions outside the polygon from the overlapping triangles with opposite orientation as one goes around the polygon following the points in the clockwise orientation. In fact this probably means that even for a general polygon which is not convex, including highly irregular shapes, the same correcting cancellation should work, so that the same procedure can be used in that case, for any chosen reference spoke point for the triangulation. We should check this guess before relying on it, no?the geometric center of philly?Klaus, Would you like me to incorporate your worksheet here and we coauthor this worksheet?