In this article
1. Theory
The centroid, \(\bar{r}\) , of a solid is defined as, $$ \bar{r}=\frac{1}{\text{Volume}}\iiint_{V}r\ dV=\frac{1}{\text{Volume}}\iiint_{V}\left(x,y,z\right)\ dV $$
where, \(r\) is a position vector. The centroid has a few very interesting properties such as:
⤷(1)
- For a homogeneous solid the center of mass is located at the centroid
- For any linear function, the average of the function inside the solid can be calculated by simply evaluating the function at the centroid point.
$$ \bar{r}=\frac{1}{\text{Volume}}\iiint_{V}\left(x,y,z\right)\ dV=\frac{1}{\text{Volume}}\iiint_{V}\left(x\frac{\partial x}{\partial x},y\frac{\partial y}{\partial y},z\frac{\partial z}{\partial z}\right)\ dV $$
⤷(2)
$$ \bar{r}=\frac{1}{\text{Volume}}\frac{1}{2}\iiint_{V}\left(\frac{\partial x^{2}}{\partial x},\frac{\partial y^{2}}{\partial y},\frac{\partial z^{2}}{\partial z}\right)\ dV $$
⤷(3)
$$ \bar{x}=\frac{1}{\text{Volume}}\frac{1}{2}\iint_{A}x^{2}\ n_{x}\ dA $$
⤷(4)
$$ \bar{y}=\frac{1}{\text{Volume}}\frac{1}{2}\iint_{A}y^{2}\ n_{y}\ dA $$
⤷(5)
$$ \bar{z}=\frac{1}{\text{Volume}}\frac{1}{2}\iint_{A}z^{2}\ n_{z}\ dA $$
where, \(\bar{x}\) , \(\bar{y}\) and \(\bar{z}\) are the components of the centroid in \(x\) , \(y\) and \(z\) directions respectively. The unit normal to the enclosing surface, \(\hat{n}=\left(n_{x},n_{y},n_{z}\right)\) , points outward of the solid. Thus, we can evaluate the centroid by calculation of second moment of area over the enclosing surface of the solid.
⤷(6)
If we divide the enclosing surface into triangles, then for each of the triangle the normal remains constant and therefore can be moved out of the integral. Let us consider the integration of \(x\text{-coordinate}\) for an individual triangle in space.
$$ I_{x}=\iint_{A}x^{2}\ n_{x}\ dA $$
Since \(n_{x}\) is constant the integration becomes,
⤷(7)
$$ I_{x}=n_{x}\iint_{A}x^{2}\ dA $$
The area integration needs to be carried out in the plane of the triangular element. Let us consider an arbitrary triangle with vertex 0, 1, 2. Let us also consider a new 2-dimensional \(u-v\) coordinate system with origin at the vertex 0 of the triangle, the \(u\text{-axis}\) aligned to base of the triangle pointing towards vertex 1. The \(v\text{-axis}\) is perpendicular to the \(u\text{-axis}\) in the plane of the triangle. The integration therefore becomes,
⤷(8)
$$ I_{x}=n_{x}\iint_{A}x^{2}\ du\, dv $$
The integration can be vastly simplified by using natural coordinate system, \(\xi-\eta\) , such that the \(\xi\text{-axis}\) is aligned to \(u\text{-axis}\) and scaled such that the vertex 0 is at \(\xi=0\) and vertex 1 is at \(\xi=1\) . The \(\eta\text{-axis}\) is aligned to triangle edge 0-2 and scaled such that the vertex 0 is at \(\eta=0\) and vertex 2 is at \(\eta=1\) . The region of integration in the \(\xi-\eta\) coordinate system therefore simplifies to a right-angled triangle, which is easy to integrate. The integrand, \(x^{2}\ du\, dv\) , needs to be transformed for the \(\xi-\eta\) coordinate system. The linear transformation for \(x\) can be obtained as,
⤷(9)
$$ x=c_{0}+c_{1}\xi+c_{2}\eta $$
Now, substituting the constraints on the coordinate system, at \(\xi=0,\eta=0\) ; \(x=x_{0}\) , we get,
⤷(10)
$$ x_{0}=c_{0}+c_{1}\times0+c_{2}\times0 $$
⤷(11)
$$ c_{0}=x_{0} $$
Substituting, at \(\xi=1,\eta=0\) ; \(x=x_{1}\) , we get,
⤷(12)
$$ x_{1}=x_{0}+c_{1}\times1+c_{2}\times0 $$
⤷(13)
$$ c_{1}=x_{1}-x_{0} $$
Substituting, at \(\xi=0,\eta=1\) ; \(x=x_{2}\) , we get,
⤷(14)
$$ x_{2}=x_{0}+c_{1}\times0+c_{2}\times1 $$
⤷(15)
$$ c_{2}=x_{2}-x_{0} $$
Therefore in the new coordinate system,
⤷(16)
$$ x=x_{0}+\left(x_{1}-x_{0}\right)\xi+\left(x_{2}-x_{0}\right)\eta $$
similarly,
⤷(17)
$$ y=y_{0}+\left(y_{1}-y_{0}\right)\xi+\left(y_{2}-y_{0}\right)\eta $$
and,
⤷(18)
$$ z=z_{0}+\left(z_{1}-z_{0}\right)\xi+\left(z_{2}-z_{0}\right)\eta $$
The integration now becomes,
⤷(19)
$$ I_{x}=n_{x}\iint_{A}\left[x_{0}+\left(x_{1}-x_{0}\right)\xi+\left(x_{2}-x_{0}\right)\eta\right]^{2}\ du\, dv $$
Also elemental area,
⤷(20)
$$ du\, dv=\frac{1}{\left|\partial\left(\xi,\eta\right)/\partial\left(u,v\right)\right|}d\xi\, d\eta $$
where, \(\left|\partial\left(\xi,\eta\right)/\partial\left(u,v\right)\right|\) is the determinant of the Jacobian. In this linear transformation it is just the ratio of area of the triangle in the \(\xi-\eta\) coordinate system to area of triangle in \(u-v\) coordinate system. Since the area of the triangle in \(\xi-\eta\) coordinate system is 0.5, the scaling factor becomes,
⤷(21)
$$ \left|\vec{A}_{T}\right|=\frac{1}{\left|\partial\left(\xi,\eta\right)/\partial\left(u,v\right)\right|}=\left|\overrightarrow{0-1}\times\overrightarrow{0-2}\right| $$
where \(\vec{A}_{T}=\overrightarrow{0-1}\times\overrightarrow{0-2}\) is the cross-product of the vectors made from edges 0-1 and 0-2. Therefore \(\left|\vec{A}_{T}\right|\) is twice the area of the triangle. The integration therefore becomes,
⤷(22)
$$ I_{x}=\left|\vec{A}_{T}\right|\, n_{x}\iint_{A}\left[x_{0}+\left(x_{1}-x_{0}\right)\xi+\left(x_{2}-x_{0}\right)\eta\right]^{2}\ d\xi\, d\eta $$
with limits of the integration,
⤷(23)
$$ I_{x}=\left|\vec{A}_{T}\right|\, n_{x}\intop_{\eta=0}^{\eta=1}\ \intop_{\xi=0}^{\xi=1-\eta}\left[x_{0}+\left(x_{1}-x_{0}\right)\xi+\left(x_{2}-x_{0}\right)\eta\right]^{2}\ d\xi\ d\eta $$
Using Maxima symbolic manipulator we can simplify this as,
⤷(24)
$$ I_{x}=\left|\vec{A}_{T}\right|\, n_{x}\left(\frac{x_{0}^{2}+x_{1}^{2}+x_{2}^{2}+x_{0}x_{1}+x_{0}x_{2}+x_{1}x_{2}}{12}\right) $$
Therefore the formula for centroid of a solid with \(N\) enclosing triangles becomes,
⤷(25)
$$ \bar{x}=\frac{1}{\text{Volume}}\frac{1}{2}\sum_{i=0}^{N-1}\left|\vec{A}_{Ti}\right|\, n_{xi}\left(\frac{x_{0}^{2}+x_{1}^{2}+x_{2}^{2}+x_{0}x_{1}+x_{0}x_{2}+x_{1}x_{2}}{12}\right)_{i} $$
⤷(26)
$$ \bar{x}=\frac{1}{24}\frac{1}{\text{Volume}}\sum_{i=0}^{N-1}\left|\vec{A}_{Ti}\right|\, n_{xi}\left(x_{0}^{2}+x_{1}^{2}+x_{2}^{2}+x_{0}x_{1}+x_{0}x_{2}+x_{1}x_{2}\right)_{i} $$
Similarly,
⤷(27)
$$ \bar{y}=\frac{1}{24}\frac{1}{\text{Volume}}\sum_{i=0}^{N-1}\left|\vec{A}_{Ti}\right|\, n_{yi}\left(y_{0}^{2}+y_{1}^{2}+y_{2}^{2}+y_{0}y_{1}+y_{0}y_{2}+y_{1}y_{2}\right)_{i} $$
and
⤷(28)
$$ \bar{z}=\frac{1}{24}\frac{1}{\text{Volume}}\sum_{i=0}^{N-1}\left|\vec{A}_{Ti}\right|\, n_{zi}\left(z_{0}^{2}+z_{1}^{2}+z_{2}^{2}+z_{0}z_{1}+z_{0}z_{2}+z_{1}z_{2}\right)_{i} $$
The above expressions for centroid eq. (27), eq. (28) and eq. (29) can be further simplified by realizing that \(\left|\vec{A}_{T}\right|\,\left(n_{x},n_{y},n_{z}\right)=\overrightarrow{0-1}\times\overrightarrow{0-2}\) , which can be calculated only once and the \(x\text{-}\) \(y\text{-}\) and \(z\text{-components}\) can be used in respective equations. It is advisable to translate the solid close to the origin for the numerical calculations so that the numerical errors due to squaring of large numbers is reduced. The actual centroid can be obtained by translating back by the same amount.
⤷(29)
2. Code
A sample program written to implement the algorithm is shown in <Code 1>
Code 1: Code for computing centroid of any 3D object
/**
* @author Sourabh Bhat (https://spbhat.in)
* Calculates the centroid point of a solid formed by a set of triangles.
* The triangles must have points either all-clockwise or
* all-anti-clockwise, looking from outside of the solid.
*
* @param triangles array of triangles
* @return centroid point of solid region enclosed by the set of triangles
*/
public class Centroid {
public static Point centroid(TriGeom [] triangles) {
Point translateBy = boundingBox(triangles).midPoint();
TriGeom [] newTriangles = translateTriangles(triangles,
-translateBy.x, -translateBy.y, -translateBy.z);
double vol = signedVolume(newTriangles);
Vector centroidPos = new Vector (0, 0, 0);
for (TriGeom t : newTriangles) {
centroidPos = centroidPos.add(scaledCentroidUnder(t));
}
centroidPos = centroidPos.mult(1.0 / vol / 24.0);
centroidPos = centroidPos.add(new Vector (
translateBy.x, translateBy.y, translateBy.z));
return new Point (centroidPos.x, centroidPos.y, centroidPos.z);
}
private static Vector scaledCentroidUnder(TriGeom tri) {
Point [] p = tri.points();
Point p0 = p[0];
Point p1 = p[1];
Point p2 = p[2];
Vector v1 = new Vector (p0, p1);
Vector v2 = new Vector (p0, p2);
Vector areaVector = v1.cross(v2);
double xc = areaVector.x * (p0.x * p0.x + p1.x * p1.x + p2.x * p2.x
+ p0.x * p1.x + p0.x * p2.x + p1.x * p2.x);
double yc = areaVector.y * (p0.y * p0.y + p1.y * p1.y + p2.y * p2.y
+ p0.y * p1.y + p0.y * p2.y + p1.y * p2.y);
double zc = areaVector.z * (p0.z * p0.z + p1.z * p1.z + p2.z * p2.z
+ p0.z * p1.z + p0.z * p2.z + p1.z * p2.z);
return new Vector (xc, yc, zc);
}
}
* @author Sourabh Bhat (https://spbhat.in)
* Calculates the centroid point of a solid formed by a set of triangles.
* The triangles must have points either all-clockwise or
* all-anti-clockwise, looking from outside of the solid.
*
* @param triangles array of triangles
* @return centroid point of solid region enclosed by the set of triangles
*/
-translateBy.x, -translateBy.y, -translateBy.z);
centroidPos = centroidPos.add(scaledCentroidUnder(t));
}
centroidPos = centroidPos.mult(1.0 / vol / 24.0);
centroidPos = centroidPos.add(
translateBy.x, translateBy.y, translateBy.z));
}
+ p0.x * p1.x + p0.x * p2.x + p1.x * p2.x);
+ p0.y * p1.y + p0.y * p2.y + p1.y * p2.y);
+ p0.z * p1.z + p0.z * p2.z + p1.z * p2.z);
}
}
The Vector and Point classes can be implemented with basic required functionality. The boundingBox() and translateTriangles() methods can be implemented to reduce error, or they can be omitted and the corresponding lines in the code can be removed. The signedVolume() method can be implemented as given in my blog on Volume calculation. For more elaborate implementations have a look at the GeometryHelper class in my CFDSolver repository.
The pdf version of the article can be downloaded using this link.