Capsule Inertia Tensor

Published January 12, 2015 by Bojan Lovrovic, posted by Bojanovski
Do you see issues with this article? Let us know.
Advertisement
Although there are robust algorithms that calculate body mass, center of mass and inertia tensor based on input triangle mesh, there are several reasons why alternatives should be used when possible. Probably the first one is speed. When working with triangles we end up with O(n) algorithm, as we need to process every of n triangles. Calculating mathematically defined body such as capsule gets us O(1) algorithm. Secondly, often we want to create a rigid body for our physics simulation but have only one mesh, a very detailed one for rendering. In such cases it is inconvenient to construct a special mesh for physics calculations and would be much easier to pass a few arguments to a function for a specific geometric primitive. Last and not so significant (at least for games) is accuracy. It is reasonable to say that a mesh that consists of triangles will never truly have round parts of its surface. This goes for spheres, cylinders, capsules and many other primitives. So if we pass our capsule as a mesh we will end up with slightly incorrect results. The analytical method that will be presented here provides accurate results. It is assumed that the reader has insight in theory behind inertia tensors and/or has an existing physics system that can make use of such a matrix. Knowledge in multivariate calculus is also assumed.

Concept

Firstly, what is meant by saying capsule? Well, simply put, it is a cylinder that (instead of flat ends) has hemispherical ends. By cylinder I mean a circular one (not elliptical). Hemispheres have the same radius. It should also be mentioned that the given body is solid (not hollow).

Capsule definition.png Figure 1. Capsule definition.

Where H is height of the cylinder part and R is radius of the hemisphere caps. This is the coordinate system that has its origin in the center of the mass of the capsule and whose axes are the axes which we will be calculating moments of inertia for.

Capsule decomposition.png Figure 2. Capsule decomposition.

The approach that we will be using to calculate the inertia tensor can be seen from Figure 2. The body will be decomposed into three parts: upper hemisphere, lower hemisphere and cylinder. Inertia tensors will be calculated for each of them and then just summed together to result in a complete capsule tensor. Here are given, in Cartesian coordinates, the well-known equations for moments and products of inertia needed for the inertia tensor in general.

eq1.gif

(eq. 1)

V represents the three dimensional region of integration and dm is an infinitesimal amount of mass at some point in our capsule. The integrand is simply squared distance from the axis.

eq2.gif

(eq. 2)

Note that in eq. 2 products of inertia are equal to zero. There exists a theorem that says that it is possible to find three (mutually perpendicular) axes for every body so the tensor produced from those axes has zero products of inertia. Those axes are called the principal axes. Furthermore, for bodies with constant density an axis of rotational symmetry is a principal axis. We have three perpendicular axes of rotational symmetry for our capsule. First one is the y axis of continuous rotational symmetry and then there are x and z axes of discrete rotational symmetry. I will not go into detail about what axes of rotational symmetry are, but I recommend to go and check it out.

eq3.gif

(eq. 3)

Since we will integrate over space eq. 3 will be used. If density is required to be some other constant than 1, the end result can simply be multiplied by that constant.

Cylinder

Before we dive into integration here is the Cartesian coordiante system we will be using for cylinder:

cylinder.png Figure 3. Cylinder and its position relative to coordinate system.

Let's start with easiest to compute moment of inertia, the one around the y axis. As previously mentioned, integrand is just a squared distance from the axis. With that in mind we can leave the Cartesian coordinate system and use the cylindrical one. The y axis remains unchanged. Also note that solutions will be written immediately without showing any steps.

eq4.gif

(eq. 4)

First we add up (integrate) over 2*pi*r interval to get the moment of inertia of a circular curve around the y axis. After that we integrate that curve over the interval R to get the moment of inertia of a disc surface perpendicular to the y axis. Lastly those discs are added up across [-H/2, H/2] interval to get the moment of inertia of a whole cylinder. The integral for the mass of a cylinder is very similar and follows the same intuition. It is given by:

eq5.gif

(eq. 5)

Now, slightly more complicated integrals are the ones for moments of inertia around x and z axis. Those are equal, so only one calculation is needed. Note that we are back in Cartesian coordinates for this one and the integral is set up like we are calculating Ixx.

eq6.gif

(eq. 6)

This is it for the cylinder as I will not go into detail about intuition for this integral since there is more work to be done for hemisphere caps.

Hemisphere caps

The following is how the hemisphere is defined relative to Cartesian coordinate system for the purpose of integration:

hemishpere.png Figure 4. Hemisphere and its position relative to coordinate system.

The easiest to calculate moment of inertia of hemispheres is also the one around y axis.

eq7.gif

(eq. 7)

We are in the cylindrical coordinate system and the intuition is similar to eq. 4. The only difference is that now the radius of a disc we need to integrate is a variable of the y position of a disc (sqrt(R^2 - y^2)). Then we just integrate those discs from 0 to R and we end up with the whole hemisphere region. As for the cylinder, the integral for the mass of a hemisphere follows the same intuition and is given by:

eq8.gif

(eq. 8)

Now, when we try to calculate the moments of inertia around x and z axis (which are equal), things get tricky. Note that integral is set up like we are calculating Ixx:

eq9.gif

(eq. 9)

hemisphere_CM.png Figure 5. Center of mass of a hemisphere.

As it can be seen in figure 5, center of mass is not in the origin of the coordinate system, but the axis for which we calculated the moment of inertia is. In order to work with this moment of inertia (to apply Steiner's rule) we need it to be in the center of mass. This is where we calculate the distance b and use it in reverse Steiner's rule:

eq10.gif

(eq. 10)

eq11.gif

(eq. 11)

Where Icm is moment of inertia of an axis that goes through the center of mass and m is hemisphere mass from eq. 8. Now it's finally possible to calculate Ixx and Izz for our hemisphere caps. We can now apply Steiner's rule to translate the axes from the center of mass of a hemisphere to the center of mass of the whole body (capsule):

eq12.gif

(eq. 12)

Moments of inertia for the hemisphere cap on the bottom of the capsule do not differ to those here.

The tensor

eq13.gif

(eq. 13)

Equations in eq. 13 give us separated masses of cylinder and hemisphere as well as the mass of the whole body. Notice how mass of the hemisphere (mhs) is multiplied by two. This is, of course, because there are two hemispheres on both ends of the capsule. Also, note our previous assumption in which density is equal to one unit of measurement, which is why the mass is equal to volume. If needed mcy and mhs can simply be multiplied by some constant density and the results will be accurate. The given code is written in such way. The full tensor uses this values and is given by:

eq14.gif

(eq. 14)

Although previously shown formulas for moments of inertia are not given by mcy and mhs, by substituting those values equality can be easily shown. The parts that represent moments of inertia of a hemisphere are multiplied by two, because we have two of them.

Code

The following code implementation is given in C-like language and is fairly optimized. #define PI 3.141592654f #define PI_TIMES2 6.283185307f const float oneDiv3 = (float)(1.0 / 3.0); const float oneDiv8 = (float)(1.0 / 8.0); const float oneDiv12 = (float)(1.0 / 12.0); void ComputeRigidBodyProperties_Capsule(float capsuleHeight, float capsuleRadius, float density, float &mass, float3 &centerOfMass, float3x3 &inertia) { float cM; // cylinder mass float hsM; // mass of hemispheres float rSq = capsuleRadius*capsuleRadius; cM = PI*capsuleHeight*rSq*density; hsM = PI_TIMES2*oneDiv3*rSq*capsuleRadius*density; // from cylinder inertia._22 = rSq*cM*0.5f; inertia._11 = inertia._33 = inertia._22*0.5f + cM*capsuleHeight*capsuleHeight*oneDiv12; // from hemispheres float temp0 = hsM*2.0f*rSq / 5.0f; inertia._22 += temp0 * 2.0f; float temp1 = capsuleHeight*0.5f; float temp2 = temp0 + hsM*(temp1*temp1 + 3.0f*oneDiv8*capsuleHeight*capsuleRadius); inertia._11 += temp2 * 2.0f; inertia._33 += temp2 * 2.0f; inertia._12 = inertia._13 = inertia._21 = inertia._23 = inertia._31 = inertia._32 = 0.0f; mass = cM + hsM * 2.0f; centerOfMass = {0.0f, 0.0f, 0.0f}; } It would seem that centerOfMass is useless, but I have left it here because it is important part of rigid body properties and often a requirement for physics engines.

Conclusion

Capsules defined in a way which was presented here are very commonly used in games and physical simulations, so this tensor should be of use to anyone building their own 3D game engine or similar software. Although, right-handed Cartesian coordinate system was used in this article, you can also use the code in left-handed system as long as y axis is the longitudinal axis (the axis that goes through caps). The somewhat tedious processes of integration were skipped and only the end results were shown. This is mostly because it would take too much space. If you are curious whether those expressions are correct, you can try them on Wolfram or try to integrate them yourself. Also, you can follow me on Twitter.

Article Update Log

December 1, 2014: Initial release December 3, 2014: Updated figures, equations and text
Cancel Save
0 Likes 7 Comments

Comments

WFP

Hi Bojan,

Nice article! I have a few questions that you can probably clear up for me.

First, it seems that the 'float3 &cm' argument is not being used for anything, and is just being hardcoded to 0 at the end of the method. Is this intentional? If so, is there a point to actually having it?

Second, am I correct in understanding that for capsules of varying masses, we would simply multiply a value representing the density by the result stored in 'float &mass'? Or is there a bit more to it than that? In other words, using the above code, how do I distinguish a representation of an object with mass 100 kg versus an object with mass 20 kg?

Thanks,

WFP

December 02, 2014 03:42 PM
Buckeye

Very nice article overall.

However, I think there is some clarification needed for technical soundness and clarity.

A. In several places, you describe integrals for the moments of inertia over in terms of Cartesian coordinates. You then switch to cylindrical coordinates for your calculations. Although it can be assumed that you intend the y-axis to be along the "vertical" axis of the cylinder through the midpoints of the caps, you should state and/or illustrate that explicitly. This comment is applicable to the hemisphere calculations also.

I.e., yes, it's easy to assume what you intend, but it takes only a few words for you to be explicit.

B. Figure 1: a minor point: it appears you intend "h" and "r" to represent the dimensions of the cylinder. However, your integrals use limits defined in terms of "H" and "R." For clarity, and to accurately tie the figure to the equations, you need to use "h" and "r," or "H" and "R," as dimensions. I suggest using "H" and "R," and clearly label those in figure 1 as dimensions, to clearly differentiate them from the axes. See comments C, D and E below, regarding your use of "H" and "R," and the assumed origin of the reference frame.

C. As mentioned in comment B above, for clarity, you need to define terms before they're used. E.g., H is defined after equation 12, although it's used in several equations prior to eq. 12.

D. Figure 1: a major point: although it appears you intend "h" and "r" as dimensions, that figure also implies that you're defining the reference frame origin as the center of the bottom cap. As your equations calculate moments of inertia about the center of the cylinder, you need to explicitly define the origin and axes with respect to the cylinder and hemispheres corresponding to those equations, with words, illustration, or both.

I.e., please make it explicit in some way that you're calculating the inertial tensor about the center of the cylinder, not the end cap.

E. Early in the article, you define R as a three dimension region, but use it as a scalar integration limit in several equations. To avoid confusion, I would suggest you use different terms in those different applications.

F. You use the term "mass" in eq. 13. I may be mistaken, but, as the units appear to indicate volume, it appears you've made an assumption about the density, a mass of 1, or a normalized mass. Can you clarify?

G. Please avoid expressions such as "obvious," "of course," or "very intuitive." Those reflect your personal opinions, and not necessarily those of the reader. :-)

H. "There exists a theorem ..." If you can, please name the theorem or provide a link to a proof. I personally know your statement to be true, but am unable to find a good reference.

I. In the code example, I would suggest you do not define constants or variables beginning with a number. E.g., perhaps change "2PI" to "PItimes2," "TwoPi" or similar. Just a programming practice.

Again, good article. Thanks for your work.

December 02, 2014 04:16 PM
Bojanovski

Hi Bojan,

Nice article! I have a few questions that you can probably clear up for me.

First, it seems that the 'float3 &cm' argument is not being used for anything, and is just being hardcoded to 0 at the end of the method. Is this intentional? If so, is there a point to actually having it?

Second, am I correct in understanding that for capsules of varying masses, we would simply multiply a value representing the density by the result stored in 'float &mass'? Or is there a bit more to it than that? In other words, using the above code, how do I distinguish a representation of an object with mass 100 kg versus an object with mass 20 kg?

Thanks,

WFP

Thanks for pointing cm (center of mass) out. The new version of the article has that part a lot clearer now.

As for your question about density, the answer is yes, it is as simple as multiplying the masses. There is nothing more to it as long as it is constant. The new version of the article has that part also explained better.

Thank you,

Bojan

December 03, 2014 01:50 AM
Bojanovski

Very nice article overall.

However, I think there is some clarification needed for technical soundness and clarity.

A. In several places, you describe integrals for the moments of inertia over in terms of Cartesian coordinates. You then switch to cylindrical coordinates for your calculations. Although it can be assumed that you intend the y-axis to be along the "vertical" axis of the cylinder through the midpoints of the caps, you should state and/or illustrate that explicitly. This comment is applicable to the hemisphere calculations also.

I.e., yes, it's easy to assume what you intend, but it takes only a few words for you to be explicit.

B. Figure 1: a minor point: it appears you intend "h" and "r" to represent the dimensions of the cylinder. However, your integrals use limits defined in terms of "H" and "R." For clarity, and to accurately tie the figure to the equations, you need to use "h" and "r," or "H" and "R," as dimensions. I suggest using "H" and "R," and clearly label those in figure 1 as dimensions, to clearly differentiate them from the axes. See comments C, D and E below, regarding your use of "H" and "R," and the assumed origin of the reference frame.

C. As mentioned in comment B above, for clarity, you need to define terms before they're used. E.g., H is defined after equation 12, although it's used in several equations prior to eq. 12.

D. Figure 1: a major point: although it appears you intend "h" and "r" as dimensions, that figure also implies that you're defining the reference frame origin as the center of the bottom cap. As your equations calculate moments of inertia about the center of the cylinder, you need to explicitly define the origin and axes with respect to the cylinder and hemispheres corresponding to those equations, with words, illustration, or both.

I.e., please make it explicit in some way that you're calculating the inertial tensor about the center of the cylinder, not the end cap.

E. Early in the article, you define R as a three dimension region, but use it as a scalar integration limit in several equations. To avoid confusion, I would suggest you use different terms in those different applications.

F. You use the term "mass" in eq. 13. I may be mistaken, but, as the units appear to indicate volume, it appears you've made an assumption about the density, a mass of 1, or a normalized mass. Can you clarify?

G. Please avoid expressions such as "obvious," "of course," or "very intuitive." Those reflect your personal opinions, and not necessarily those of the reader. :-)

H. "There exists a theorem ..." If you can, please name the theorem or provide a link to a proof. I personally know your statement to be true, but am unable to find a good reference.

I. In the code example, I would suggest you do not define constants or variables beginning with a number. E.g., perhaps change "2PI" to "PItimes2," "TwoPi" or similar. Just a programming practice.

Again, good article. Thanks for your work.

Whoa, that's a lot of work you gave me, but your notes are sound and I did my best to update the article. If there is anything else, let me know.

Thanks,

Bojan

December 03, 2014 01:54 AM
Buckeye

Good job, Bojan. Thank you. Thumbs up for the article.

One more comment: you need to change all occurrences in the code of the defined constant "2PI" to "PI_TIMES2." I.e., you missed:


hsM = 2PI*oneDiv3*rSq*capsuleRadius*density;
// should be
hsM = PI_TIMES2*oneDiv3*rSq*capsuleRadius*density;
December 03, 2014 03:16 AM
Dave Hunt

Nice work, Bojan!

January 12, 2015 05:38 AM
MrRowl

Thanks for this. However, eq 14 has an error! Where it says H^2/2 that should be (H/2)^2.

The code is correct, but of course nobody would just pinch that :)

October 21, 2017 11:19 PM
You must log in to join the conversation.
Don't have a GameDev.net account? Sign up!

This article tries to simplify the math behind calculating inertia tensor of a cylinder as much as possible without skipping on any important details. Implementation in C/C++ is presented at the end.

Advertisement

Other Tutorials by Bojanovski

Advertisement