# Finding the centroid of the area between two curves

When I have an area bounded by curves, is there a built-in way to find the center of the area? Or do I have to plot it first and then use ComponentMeasurements on it?

For example: the area under $y= 1-x^2/4$ and above the $x$ axis:

-
The centroid is given by $\bar x = (1/A) \iint x \; dA$, $\bar y = (1/A) \iint y \; dA$, where $A=\iint \; dA$ and the integrals are over the area. Does that answer your question, or is it about how to set up such integrals in Mathematica? – Michael E2 Jan 2 at 14:58
I know how to calculate it with integrals, I want to know if there is a function to do it directly (without me having to calculate $dA$, $X_e$, $Y_e$ first) – a3f Jan 2 at 15:14
then you´ll have to specify more precisely what your input looks like – Yves Klett Jan 2 at 15:31
@YvesKlett My input looks like I wrote above: A number of curves which enclose an area – a3f Jan 2 at 15:51
@YvesKlett I thought it was clear enough... His question basically is: "Is there a built in that takes a curve (or a set of curves) and finds the centroid of the enclosing area using calculus or should I plot it, rasterize it and use ComponentsMeasurements". Looks like they already know what Comp..M... does – rm -rf Jan 2 at 16:09
show 6 more comments

This is inelegent for the specific example, but may be useful for more greneral cases:

boolf[x_, y_] := (y > 0 && y < 1 - x^2/4);
area = NIntegrate[If[  boolf[x, y] , 1, 0], {x, -20, 20}, {y, -20, 20}]
-> 2.66667
cx = NIntegrate[If[  boolf[x, y] , x, 0], {x, -20, 20}, {y, -20, 20}]/ area
-> ~10^-16
cy = NIntegrate[If[  boolf[x, y] , y, 0], {x, -20, 20}, {y, -20, 20}]/area
-> 0.4


The exact results for the example are by the way:

area = Integrate[( 1 - x^2/4 ) , {x, -2, 2}]  -> 8/3
cy = Integrate[( 1 - x^2/4 )/2  ( 1 - x^2/4 ) , {x, -2, 2}] / area  -> 2/5
cx = Integrate[ x  ( 1 - x^2/4 ) , {x, -2, 2}] / area   -> 0


Edit: further explintion: we construct a logical function of x,y which is True inside your region and False otherwise, then perform area integrals over the entire plane using the definitions of area and centriod, and relying on the integrands to be zero outside the region of interest.

As noted by Daniel in a comment that seems to have dissapeard(?) you can use the analytic Integrate function over the entire plane (+- Infinity):

fa[x_,y_]:= If[(y > 0 &&  y < 1 - x^2/4), 1, 0];
area = Integrate[fa[x, y], {x, -Infinity, Infinity}, {y, -Infinity, Infinity}]
-> 8/3
cx = Integrate[ x fa[x,y], {x, -Infinity, Infinity}, {y, -Infinity, Infinity}]/ area
-> 0
cy = Integrate[ y fa[x,y], {x, -Infinity, Infinity}, {y, -Infinity, Infinity}]/area
-> 2/5


This does need a bit of warning, Integrate[] fails for even slightly more complicated expressions, so NIntegrate[] is a bit safer (though it yields a numerical approximation).

-
oops worked on that before i saw Danial's note. And it does indeed work with plain Integrate (which surprised me..) – george2079 Jan 2 at 17:30
Thanks, But could u explain what exactly you did ? (I'm new to Mathematica) – a3f Jan 2 at 17:55
@a3f: What is it you still don't understand? – murray Jan 2 at 19:17
The code is a perhaps more transparent if you use the built-in function Boole, which already returns values of 1 (when its argument evalutaes to True and 0 otherwise). Thus: indicator[x_, y_] := Boole[y > 0 && y < 1 - x^2/4] . And then area = Integrate[indicator[x, y], {x, -20, 20}, {y, -20, 20}], cx = Integrate[x indicator[x, y], {x, -20, 20}, {y, -20, 20}]/area, etc. – murray Jan 2 at 19:20
As you know bounds on the region, you might as well use them, i.e., integrate just for {x, -2, 2} and {y, 0, 1}. Especially in case you integrate numerically. – murray Jan 2 at 19:42
show 2 more comments

I finally found the following, which I had done for someone who had wanted to find centroids of polygons, for which NIntegrate seems too slow. I thought it was pretty cool at the time and wanted to share it. I guess it's pretty clear that the simple answer to the original question of @a3f is that there is no built-in centroid function. I hope you will agree that this is relevant enough to the question to be worth posting.

The math behind the following is based the ideas in @Jens hints. If you would like an explanation, ask and I will append one.

AreaIntegrand = Compile[{Px, Py, Qx, Qy},
(Qx + Px) (Qy - Py)];
XMomentIntegrand = Compile[{Px, Py, Qx, Qy},
1/6 (Px - Qx) (2 Px Py + Py Qx + Px Qy + 2 Qx Qy)];
YMomentIntegrand = Compile[{Px, Py, Qx, Qy},
-(1/6) (Py - Qy) (2 Px Py + Py Qx + Px Qy + 2 Qx Qy)];
Area = Plus @@ AreaIntegrand @@@ Join @@@ #/2 &;
XMoment = Plus @@ XMomentIntegrand @@@ Join @@@ # &;
YMoment = Plus @@ YMomentIntegrand @@@ Join @@@ # &;
XBar = XMoment[#]/Area[#] &;
YBar = YMoment[#]/Area[#] &;
Centroid[poly_Polygon] := Centroid[First[poly]];
Centroid[poly_List?(Depth[#] == 4 &)] := {#[[1]], #[[2]]}/#[[3]] &@
((#.(Sign /@ Last[#]) &)@({XMoment /@ #, YMoment /@ #, Area /@ #} &)@
(Partition[#, 2, 1, 1] & /@ poly));
Centroid[poly_List?(Depth[#] == 3 &)] :=
{XMoment[#], YMoment[#]} / Area[#] &[Partition[poly, 2, 1, 1]];


Beware: The sub-functions (Area, XMoment, XBar, etc.) assume that the argument is a list of the edges of a polygon (Partition[poly, 2, 1, 1]).

As a test I loaded Italy and found its centroid:

itPoly = CountryData["Italy", "Polygon"];
Centroid[itPoly]
(* {12.0822, 42.7896} *)


It's the centroid of a "flat-earth" Italy, since no adjustment for spherical coordinates is made.

Graphics[{itPoly, PointSize[Large], Red, Point[Centroid[itPoly]]}]


It's not too slow:

Centroid[itPoly] // AbsoluteTiming
(* {0.007083, {12.0822, 42.7896}} *)


The method is exact on polygons; if one wants exact results, Compile needs to be replaced by Function. Compile seems to speed things up by a factor of 2. Centroid is faster than NIntegrate on polygons:

{#[[1]], #[[2]]}/#[[3]] &@(Plus @@ (Module[{X, Y, a, xm, ym},
X = Interpolation[Append[#, First[#]] &[First /@ #], InterpolationOrder -> 1];
Y = Interpolation[Append[#, First[#]] &[Last /@ #], InterpolationOrder -> 1];
a = NIntegrate[X[t] Y'[t], {t, 1, 1 + Length[#]}];
xm = NIntegrate[-X[t] Y[t] X'[t], {t, 1, 1 + Length[#]}];
ym = NIntegrate[X[t] Y[t] Y'[t], {t, 1, 1 + Length[#]}];
{xm, ym, a}
] & /@ First[itPoly])) // AbsoluteTiming
(* {4.529647, {12.0822, 42.7896}} *)


Curves

The function Centroid above is not exact on curves, but it may be accurate enough for some work. In most cases, I'd bet that NIntegrate would beat it. First you have to convert the curve to an approximate polygon and then apply Centroid.

Graphs

For the example proposed by @a3f, we can extract the plotted points of the graph with Cases. Since the base is straight, they form the polygon and we can apply Centroid:

(paraPoly = Cases[Plot[1 - x^2/4, {x, -2, 2}, PlotPoints -> 800,
EvaluationMonitor :> Sow[{x, 1 - x^2/4}]],
Line[pts_] :> pts, \[Infinity]];
Centroid[paraPoly]) // AbsoluteTiming
(* {0.035418, {5.32535*10^-10, 0.4}} *)


The relative error of the y-coordinate is about $10^{-6}$.

Graphics[{Polygon[paraPoly], PointSize[Large], Red,
Point[Centroid[paraPoly]]}, Frame -> True]


If we compare with Integrate and NIntegrate as @murray, @george2079 have answered, we get:

(a = Integrate[(1 - x^2/4), {x, -2, 2}];
xm = Integrate[x (1 - x^2/4), {x, -2, 2}];
ym = Integrate[(1 - x^2/4)/2 (1 - x^2/4), {x, -2, 2}];
{xm, ym} / a) // AbsoluteTiming
(* {0.020960, {0, 2/5}} *)

(a = NIntegrate[(1 - x^2/4), {x, -2, 2}];
xm = NIntegrate[x (1 - x^2/4), {x, -2, 2}];
xm = NIntegrate[(1 - x^2/4)/2 (1 - x^2/4), {x, -2, 2}];
{xm, ym} / a) // AbsoluteTiming
(* {0.019479, {0., 0.4}} *)


The first is exact and the second is accurate to $MachinePrecision. Parametric Curves I made up, at random, simple closed curve (picture down below), f[t_] := (2 + Sin[8 \[Pi] t] + Cos[2 \[Pi] t]) {Cos[2 \[Pi] t], Sin[2 \[Pi] t]};  and found that for 11000 plot points (no recursion), the Centroid calculation took about as long as the NIntegrate. So we can compare the results. (polarPoly = Cases[ParametricPlot[Evaluate[f[t]], {t, 0, 1}, PlotPoints -> 11000, MaxRecursion -> 0], Line[l_] :> l, \[Infinity]]; Centroid[polarPoly]) // AbsoluteTiming (* {0.205172, {0.95, 6.16864*10^-12}} *) Module[{a, xm, ym}, a = NIntegrate[Evaluate[f[t][[1]] f'[t][[2]]], {t, 0, 1}]; xm = NIntegrate[Evaluate[-f[t][[1]] f[t][[2]] f'[t][[1]]], {t, 0, 1}]; ym = NIntegrate[Evaluate[f[t][[1]] f[t][[2]] f'[t][[2]]], {t, 0, 1}]; {xm, ym}/a ] // AbsoluteTiming (* {0.207635, {0.95, -1.65599*10^-16}} *) Module[{a, xm, ym}, a = Integrate[Evaluate[f[t][[1]] f'[t][[2]]], {t, 0, 1}]; xm = Integrate[-Evaluate[f[t][[1]] f[t][[2]] f'[t][[1]]], {t, 0, 1}]; ym = Integrate[ Evaluate[f[t][[1]] f[t][[2]] f'[t][[2]]], {t, 0, 1}]; {xm, ym}/a ] // AbsoluteTiming (* {9.105406, {19/20, 0}} *)  The x-coordinate of Centroid is off by less than$10^{-7}\$, which is not bad. With a complicated integrand Integrate bogs down and NIntegrate is clearly superior (surprise, surprise).

-
 Pretty slick. Easy +1. (Even though you answered quite a while ago.) – kale May 10 at 18:11

This is a nice application of Stokes' theorem. In fact it could almost be a homework problem - I hope it isn't, so I'll give only the first steps:

Observe the following results:

Curl[{0, x^2/2, 0}, {x, y, z}]

(* ==> {0, 0, x} *)

Curl[{-y^2/2, 0, 0}, {x, y, z}]

(* ==> {0, 0, y} *)

Curl[{-y/2, x/2, 0}, {x, y, z}]

(* ==> {0, 0, 1} *)


Now think about the surface integrals involved in the calculation: one is simply the area, i.e. its integrand is 1. Then there are also the integrands x and y for the Cartesian components of the centroid.

But since you're given the curves bounding the area, you would be better off writing the area integrations as line integrals along the bounding curves. This is where Stokes' theorem comes in, combined with the above results.

The rest should be left for you to complete. This method, once you've written it down, will have the added advantage that it can give you the centroid of the projection of any three-dimensional curve onto the xy plane.

-
 It's not exactly homework. I am solving some problems in the book and wondered how I could use Mathematica to check if I was right. Thanks for your answer :) – a3f Jan 2 at 19:56 The nice thing about using this theorem is that you end up having to do only half as many integrals as in the accepted answer... so it's worth reading about Stokes' theorem, too! – Jens Jan 2 at 20:38 @a3f: Presumably, then, you're just studying centroids as an application of single-variable calculus; in that case, using curl and Stokes's Theorem is something you might want to save for the future, when you study multivariable calculus. – murray Jan 2 at 20:40

Here's a more down-to-earth solution in Mathematica that's probably closer to your current level of mathematics study. Just directly use the usual one-dimensional integrals involved in finding plane area and centroids.

f[x_] := 1 - x^2/4
area = Integrate[f[x], {x, -2, 2}]
yMoment = Integrate[x f[x], {x, -2, 2}]
xMoment = Integrate[(1/2) f[x]^2, {x, -2, 2}]
centroid = {yMoment/area, xMoment/area}


Of course there's no real need to calculate the x-coordinate, since by symmetry it must be 0.

-
 We already started multi-variable calculus but thanks for your example. – a3f Jan 2 at 22:35

I had previously used the following function as part of designing this site's logo, but I suppose it would be useful to explicitly have it as an answer to this question:

PolygonCentroid[pts_?MatrixQ] := With[{dif = Map[Det, Partition[pts, 2, 1, {1, 1}]]},
ListConvolve[{{1, 1}}, Transpose[pts], {-1, -1}].dif/(3 Total[dif])]


A test, using OP's example:

reg = First[Cases[Plot[1 - x^2/4, {x, -2, 2}, PlotPoints -> 95], Line[l_] :> l, ∞]];

PolygonCentroid[reg]
{1.00974*10^-6, 0.399977}


Pretty close for government work, I think. I exploited Plot[]'s capabilities for adaptive sampling in this case to help compute a more accurate estimate for the centroid. You can tweak PlotPoints -> 95 or other options as needed for refinement.

-