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

Edit: Preface

After being on this site for a few months, I've learned quite a bit about what Mathematica can do quickly. Below I've revised my code to take advantage of built-in vectorized functions. One could compile the functions as well, but the advantages perhaps outweigh the disadvantages. The increase in speed is slight on approximate real numbers (because the number crunching is done in a very few function calls and so is already "mostly compiled" anyway); however, the first function centroid1 is exact on polygons, and a compiled function would not return an exact result on exact input. I added a new one that implements a Simpson's Rule type integration that is exact on figures whose boundary can be parametrized by quadratic polynomials. I also updated the timings, since I have a new machine.

Revised original method

Originally I had done this for someone who had wanted to find centroids of polygons, for which NIntegrate seemed 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.

The code below replaces the previous code but is based on the same mathematical ideas as before. The line integrals around the boundary, $$A={1\over2} \int x\;dy - y\;dx, \quad M_x=-\int x\,y\;dx, \quad M_y=\int x\,y\;dy\,,$ which by Green's Theorem represent the area, x moment, and y moment respectively, are approximated using a linear interpolation of the boundary at sample points.

area1[pts_] := Module[{x, y}, {x=#[[1]],y=#[[2]]} &@ Transpose @ pts;
   Total[x RotateLeft[y] - RotateLeft[x] y]] / 2;
xMoment1[pts_] := Module[{x, y}, {x=#[[1]],y=#[[2]]} &@ Transpose @ pts; 
   Total[(x - RotateLeft[x]) (2 x y + y RotateLeft[x] + 
        x RotateLeft[y] + 2 RotateLeft[x] RotateLeft[y])]] / 6;
yMoment1[pts_] := Module[{x, y}, {x=#[[1]],y=#[[2]]} &@ Transpose @ pts; 
   Total[(RotateLeft[y] - y) (2 x y + y RotateLeft[x] + 
        x RotateLeft[y] + 2 RotateLeft[x] RotateLeft[y])]] / 6;

centroid1[poly_Polygon] := centroid1[First[poly]];
centroid1[poly_List?(Depth[#] == 4 &)] := (* multiple-polygon figure (e.g. Italy) *)
   {#[[1]], #[[2]]} / #[[3]] &[ 
        # . Sign@Last[#] & @ {xMoment1 /@ #, yMoment1 /@ #, area1 /@ #} & @ poly];
centroid1[poly_List?(Depth[#] == 3 &)] := {xMoment1[#], yMoment1[#]} / area1[#] & @ poly;

Polygons

As a test I loaded Italy and found its centroid:

itPoly = CountryData["Italy", "Polygon"];
centroid1[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[centroid1[itPoly]]}]

The new code is faster:

centroid1[itPoly] // AbsoluteTiming (* new *)
(* {0.000468, {12.0822, 42.7896}} *)

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

The function centroid1 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
(* {2.616189, {12.0822, 42.7896}} *)

Curves

The function centroid1 above is not exact on curves, but it may be accurate enough for some work. In most cases, NIntegrate is hard to beat it for accuracy. To apply centroid1, first you have to convert the curve to an approximate polygon (pick sample points).

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 centroid1:

(paraPoly = Transpose@{#, 1 - #^2/4} &@Range[-2., 2, 4/1000];
 centroid1[paraPoly]) // AbsoluteTiming
(* {0.000356, {6.83717*10^-17, 0.4}} *)

The error in the y-coordinate is less than $10^{-6}$.

Graphics[{Polygon[paraPoly], PointSize[Large], Red, 
  Point[centroid1[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.020392, {0, 2/5}} *)

(a = NIntegrate[(1 - x^2/4), {x, -2, 2}];
  xm = NIntegrate[x (1 - x^2/4), {x, -2, 2}];
  ym = NIntegrate[(1 - x^2/4)/2 (1 - x^2/4), {x, -2, 2}];
  {xm, ym} / a) // AbsoluteTiming
(* {0.013740, {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_] := {Cos[2 \[Pi] t] (2 + Cos[2 \[Pi] t] + Sin[8 \[Pi] t]), Sin[2 \[Pi] t] (2 + Cos[2 \[Pi] t] + Sin[8 \[Pi] t])};

and found that for 200000 points, the centroid1 calculation took about as long as the NIntegrate. So we can compare the results.

(polarPoly = Transpose @ f[Range[0., 1, 1/200000]];
   centroid1[polarPoly]) // AbsoluteTiming
(* {0.104392, {0.95, -1.38784*10^-14}} *)

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.108268, {0.95, -1.11203*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
(* {7.968150, {19/20, 0}} *)

The errors in the x coordinate are a little over 10^-10 and 10^-16 (or $MachinePrecision) for centroid1 and NIntegrate respectively. (But on 5000 points centroid1 is 100 times faster and has an error in x of a little over 10^-7, which is good enough for plotting.)

---

Addition -- Simpson's Rule

If we have an even number of points, we can use a method similar to Simpson's Rule. Think of the points as alternating between "endpoints" and "midpoints," and then use quadratic interpolation to approximate the boundary of the region. This tends to be better on curved boundaries and worse on straight boundaries such as polygons. With polygons, we will get an inaccurate result unless we bisect the segments; so from the points of view of both accuracy and speed, it loses over the linear method of centroid1.

With[{
    areaFn = Total[(-3 #1 - 4 #2 + #3) #4 + 4 (#1 - #3) #5 - (#1 - 4 #2 - 3 #3) #6]/6 &,
    momentFn = Total[(-10 #1 - 6 #2 + #3) #4^2 + (6 #1 - 8 #2 + 2 #3) (#4 #5) + 
       8 (#1 - #3) #5^2 + (-#1 + #3) (#4 #6) + (-2 #1 + 8 #2 - 6 #3) (#5 #6) +
       (-#1 + 6 #2 + 10 #3) #6^2]/30 &},
 area2[pts_] := Module[{x, y}, {x=#[[1]],y=#[[2]]} &@ Transpose @ pts;
   areaFn[x[[1 ;; -1 ;; 2]], x[[2 ;; -1 ;; 2]], 
     RotateLeft[x[[1 ;; -1 ;; 2]]], y[[1 ;; -1 ;; 2]], 
     y[[2 ;; -1 ;; 2]], RotateLeft[y[[1 ;; -1 ;; 2]]]]
   ];
 yMoment2[pts_] := Module[{x, y}, {x=#[[1]],y=#[[2]]} &@ Transpose @ pts;
   momentFn[x[[1 ;; -1 ;; 2]], x[[2 ;; -1 ;; 2]], 
     RotateLeft[x[[1 ;; -1 ;; 2]]], y[[1 ;; -1 ;; 2]], 
     y[[2 ;; -1 ;; 2]], RotateLeft[y[[1 ;; -1 ;; 2]]]]
   ];
 xMoment2[pts_] := Module[{x, y}, {x=#[[1]],y=#[[2]]} &@ Transpose @ pts;
   -momentFn[y[[1 ;; -1 ;; 2]], y[[2 ;; -1 ;; 2]], 
      RotateLeft[y[[1 ;; -1 ;; 2]]], x[[1 ;; -1 ;; 2]], 
      x[[2 ;; -1 ;; 2]], RotateLeft[x[[1 ;; -1 ;; 2]]]]
   ]
 ];

centroid2[poly_Polygon] := centroid2[First[poly]];
centroid2[poly_List?(Depth[#] == 4 &)] /; And @@ (EvenQ /@ Length /@ poly) :=
  {#[[1]], #[[2]]} / #[[3]] &[#.(Sign@Last[#]) &[
      {xMoment2 /@ #, yMoment2 /@ #, area2 /@ #} &@ poly]];
centroid2[poly_List?(Depth[#] == 3 &)] /; EvenQ @ Length[poly] :=
  {xMoment2[#], yMoment2[#]} / area2[#] &@ poly;

It's exact on regions with boundaries that have quadratic parametrizations. Here care is taken to have the bottom segment bisected by the origin as midpoint in order to get the exact answer.

(paraPoly = Transpose@{#, 1 - #^2/4} &@ Range[-2, 2, 4/4] ~Join~ {{0, 0}};
 centroid2[paraPoly]) // AbsoluteTiming
(* {0.000289, {0, 2/5}} *)

On the parametric curve, with 1000 sample points, the error in x is a little over 10^-7, the same as centroid1 on 5000 points.

f[t_] := {Cos[2 \[Pi] t] (2 + Cos[2 \[Pi] t] + Sin[8 \[Pi] t]), 
          Sin[2 \[Pi] t] (2 + Cos[2 \[Pi] t] + Sin[8 \[Pi] t])};
(polarPoly = Transpose @ f[Range[0., 1, 1/999]];
 centroid2[polarPoly]) // AbsoluteTiming
(* {0.000714, {0.95, 5.44169*10^-10}} *)

Answer 2

Now in V10 we have ImplicitRegion and RegionCentroid to do this easily:

reg = ImplicitRegion[y > 0 && y <= 1 - x^2/4, {x, y}];

Then:

RegionCentroid[reg]

{0, 2/5}

Answer 3

This is inelegant for the specific example, but may be useful for more general 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 explanation: 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 centroid, and relying on the integrands to be zero outside the region of interest.

As noted by Daniel in a comment that seems to have disappeared(?) you can use the analytic Integrate function over the entire plane (±∞):

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

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.

Answer 5

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 6

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.