# 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 '13 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 '13 at 15:14
then you´ll have to specify more precisely what your input looks like –  Yves Klett Jan 2 '13 at 15:31
@YvesKlett My input looks like I wrote above: A number of curves which enclose an area –  a3f Jan 2 '13 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 –  The Toad Jan 2 '13 at 16:09

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 '13 at 17:30
Thanks, But could u explain what exactly you did ? (I'm new to Mathematica) –  a3f Jan 2 '13 at 17:55
@a3f: What is it you still don't understand? –  murray Jan 2 '13 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 '13 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 '13 at 19:42

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

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

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}} *)

-
Pretty slick. Easy +1. (Even though you answered quite a while ago.) –  kale May 10 '13 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 '13 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 '13 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 '13 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 '13 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.

-

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}

-