Finding the centroid of the area between two curves

Question

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:

Answer 1

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).

Answer 2

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.

Answer 3

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).

Answer 4

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.

Answer 5

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.