Calculating the Moment Of Inertia for a concave 2D polygon relative to its orgin

Question

I want to compute the moment of inertia of a (2D) concave polygon. I found this on the internet. But I'm not very sure how to interpret the formula...

Formula http://img101.imageshack.us/img101/8141/92175941c14cadeeb956d8f.gif

1) Is this formula correct?
2) If so, is my convertion to C++ correct?

float sum (0);
for (int i = 0; i < N; i++) // N = number of vertices
{
    int j = (i + 1) % N;
    sum += (p[j].y - p[i].y) * (p[j].x + p[i].x) * (pow(p[j].x, 2) + pow(p[i].x, 2)) - (p[j].x - p[i].x) * (p[j].y + p[i].y) * (pow(p[j].y, 2) + pow(p[i].y, 2));
}
float inertia = (1.f / 12.f * sum) * density;

Martijn

Looking at it again, the parts you say about cross product don't really make sense to me. [Cross product](http://en.wikipedia.org/wiki/Cross_product) is only defined for 3-D vectors, and you say you are using 2-D ones... — Larry Wang, Jul 25 '10 at 14:00
@Larry Wang: I havn't given the above formula much (any) thought, but because R^2 is isomorphic to the plane z = 0 in R^3, you can actually use cross products for 2D vectors as well, implicitly setting their z coordinate equal to zero. — Andreas Rejbrand, Jul 25 '10 at 14:06
@Andreas: Is that the solution for that problem? I encountered that problem while converting Larry's code.... — Martijn Courteaux, Jul 25 '10 at 14:12
@Martijn: As I said, I havn't studied the formula particulary carefully, but as Larry Wang pointed out, the cross product can only be computed in three dimensions. But it is common practice to use it in the plane as well, pretending the plane to be the subspace z = 0 of the three-dimensional space, i.e. adding a third z coordinate equal to zero to all planar vectors. — Andreas Rejbrand, Jul 25 '10 at 14:20
Some of the comments in your code are incorrect. The cross product always creates a vector in both 2D and 3D space. It just so happens that in 2D it's a vector that is normal to the 2D plane. It is not a scalar by any means. And there's no such thing as the cross-product of a scalar and a vector. You can multiply a vector by a scalar to get another vector, but it's not a cross-product. — duffymo, Jul 25 '10 at 14:38
@Andreas, @Martijn: Since the only reason we care about the cross product here is because we want to find its length, that seems like it would be ok to me. I don't have any idea why you're crossing a vector by a scalar though. @duffymo: It's not the scalar product. — Larry Wang, Jul 25 '10 at 14:38
The so-called "cross product with a scalar" is with the vector `xyz = (0, 0, s)`. Which matches the erroneous labeling of the vector cross-product results as a 'scalar'. — Ben Voigt, Jul 25 '10 at 15:02
@Ben: You're right, that fits perfectly. Very strange notation. — Larry Wang, Jul 25 '10 at 15:16

Larry Wang · Accepted Answer · 2010-07-25T15:15:29.830

5

#include <math.h> //for abs
float dot (vec a, vec b) {
   return (a.x*b.x + a.y*b.y);
}
float lengthcross (vec a, vec b) {
   return (abs(a.x*b.y - a.y*b.x));
}
...
do stuff
...
float sum1=0;
float sum2=0;
for (int n=0;n<N;++n)  { //equivalent of the Σ
   sum1 += lengthcross(P[n+1],P[n])* 
           (dot(P[n+1],P[n+1]) + dot(P[n+1],P[n]) + dot(P[n],P[n]));
   sum2 += lengthcross(P[n+1],P[n]);
}
return (m/6*sum1/sum2);

Edit: Lots of small math changes

edited Jul 25 '10 at 15:15

answered Jul 25 '10 at 13:51

Larry Wang

986
6
17

That's not actually the length (magnitude or two-norm) of the cross-product, that's the z-coordinate of the cross-product. Length of any vector is non-negative. Not sure which is correct, but the formula _appears_ to call for the length. – Ben Voigt Jul 25 '10 at 14:58
And, the function `length` is named wrong (it is the length squared). – Ben Voigt Jul 25 '10 at 14:59

duffymo · Answer 2 · 2010-07-25T15:27:39.270

I think you have more work to do that merely translating formulas into code. You need to understand exactly what this formula means.

When you have a 2D polygon, you have three moments of inertia you can calculate relative to a given coordinate system: moment about x, moment about y, and polar moment of inertia. There's a parallel axis theorem that allows you to translate from one coordinate system to another.

Do you know precisely which moment and coordinate system this formula applies to?

Here's some code that might help you, along with a JUnit test to prove that it works:

import java.awt.geom.Point2D;

/**
 * PolygonInertiaCalculator
 * User: Michael
 * Date: Jul 25, 2010
 * Time: 9:51:47 AM
 */
public class PolygonInertiaCalculator
{
    private static final int MIN_POINTS = 2;

    public static double dot(Point2D u, Point2D v)
    {
        return u.getX()*v.getX() + u.getY()*v.getY();
    }

    public static double cross(Point2D u, Point2D v)
    {
        return u.getX()*v.getY() - u.getY()*v.getX();
    }

    /**
     * Calculate moment of inertia about x-axis
     * @param poly of 2D points defining a closed polygon
     * @return moment of inertia about x-axis
     */
    public static double ix(Point2D [] poly)
    {
        double ix = 0.0;

        if ((poly != null) && (poly.length > MIN_POINTS))
        {
            double sum = 0.0;
            for (int n = 0; n < (poly.length-1); ++n)
            {
                double twiceArea = poly[n].getX()*poly[n+1].getY() - poly[n+1].getX()*poly[n].getY();
                sum += (poly[n].getY()*poly[n].getY() + poly[n].getY()*poly[n+1].getY() + poly[n+1].getY()*poly[n+1].getY())*twiceArea;
            }

            ix = sum/12.0;
        }

        return ix;
    }

    /**
     * Calculate moment of inertia about y-axis
     * @param poly of 2D points defining a closed polygon
     * @return moment of inertia about y-axis
     * @link http://en.wikipedia.org/wiki/Second_moment_of_area
     */
    public static double iy(Point2D [] poly)
    {
        double iy = 0.0;

        if ((poly != null) && (poly.length > MIN_POINTS))
        {
            double sum = 0.0;
            for (int n = 0; n < (poly.length-1); ++n)
            {
                double twiceArea = poly[n].getX()*poly[n+1].getY() - poly[n+1].getX()*poly[n].getY();
                sum += (poly[n].getX()*poly[n].getX() + poly[n].getX()*poly[n+1].getX() + poly[n+1].getX()*poly[n+1].getX())*twiceArea;
            }

            iy = sum/12.0;
        }

        return iy;
    }

    /**
     * Calculate polar moment of inertia xy
     * @param poly of 2D points defining a closed polygon
     * @return polar moment of inertia xy
     * @link http://en.wikipedia.org/wiki/Second_moment_of_area
     */
    public static double ixy(Point2D [] poly)
    {
        double ixy = 0.0;

        if ((poly != null) && (poly.length > MIN_POINTS))
        {
            double sum = 0.0;
            for (int n = 0; n < (poly.length-1); ++n)
            {
                double twiceArea = poly[n].getX()*poly[n+1].getY() - poly[n+1].getX()*poly[n].getY();
                sum += (poly[n].getX()*poly[n+1].getY() + 2.0*poly[n].getX()*poly[n].getY() + 2.0*poly[n+1].getX()*poly[n+1].getY() + poly[n+1].getX()*poly[n].getY())*twiceArea;
            }

            ixy = sum/24.0;
        }

        return ixy;
    }

    /**
     * Calculate the moment of inertia of a 2D concave polygon
     * @param poly array of 2D points defining the perimeter of the polygon
     * @return moment of inertia
     * @link http://www.physicsforums.com/showthread.php?t=43071
     * @link http://www.physicsforums.com/showthread.php?t=25293
     * @link http://stackoverflow.com/questions/3329383/help-me-with-converting-latex-formula-to-code
     */
    public static double inertia(Point2D[] poly)
    {
        double inertia = 0.0;

        if ((poly != null) && (poly.length > MIN_POINTS))
        {
            double numer = 0.0;
            double denom = 0.0;
            double scale;
            double mag;
            for (int n = 0; n < (poly.length-1); ++n)
            {
                scale = dot(poly[n + 1], poly[n + 1]) + dot(poly[n + 1], poly[n]) + dot(poly[n], poly[n]);
                mag = Math.sqrt(cross(poly[n], poly[n+1]));
                numer += mag * scale;
                denom += mag;
            }

            inertia = numer / denom / 6.0;
        }

        return inertia;
    }
}

Here's the JUnit test to accompany it:

import org.junit.Test;

import java.awt.geom.Point2D;

import static org.junit.Assert.assertEquals;

/**
 * PolygonInertiaCalculatorTest
 * User: Michael
 * Date: Jul 25, 2010
 * Time: 10:16:04 AM
 */
public class PolygonInertiaCalculatorTest
{
    @Test
    public void testTriangle()
    {
        Point2D[] poly =
        {
            new Point2D.Double(0.0, 0.0),
            new Point2D.Double(1.0, 0.0),
            new Point2D.Double(0.0, 1.0)
        };


        // Moment of inertia about the y1 axis
        // http://www.efunda.com/math/areas/triangle.cfm
        double expected = 1.0/3.0;
        double actual = PolygonInertiaCalculator.inertia(poly);

        assertEquals(expected, actual, 1.0e-6);
    }

    @Test
    public void testSquare()
    {
        Point2D[] poly =
        {
            new Point2D.Double(0.0, 0.0),
            new Point2D.Double(1.0, 0.0),
            new Point2D.Double(1.0, 1.0),
            new Point2D.Double(0.0, 1.0)
        };

        // Polar moment of inertia about z axis
        // http://www.efunda.com/math/areas/Rectangle.cfm
        double expected = 2.0/3.0;
        double actual = PolygonInertiaCalculator.inertia(poly);

        assertEquals(expected, actual, 1.0e-6);
    }

    @Test
    public void testRectangle()
    {
        // This gives the moment of inertia about the y axis for a coordinate system
        // through the centroid of the rectangle
        Point2D[] poly =
        {
            new Point2D.Double(0.0, 0.0),
            new Point2D.Double(4.0, 0.0),
            new Point2D.Double(4.0, 1.0),
            new Point2D.Double(0.0, 1.0)
        };

        double expected = 5.0 + 2.0/3.0;
        double actual = PolygonInertiaCalculator.inertia(poly);

        assertEquals(expected, actual, 1.0e-6);

        double ix = PolygonInertiaCalculator.ix(poly);
        double iy = PolygonInertiaCalculator.iy(poly);
        double ixy = PolygonInertiaCalculator.ixy(poly);

        assertEquals(ix, (1.0 + 1.0/3.0), 1.0e-6);
        assertEquals(iy, (21.0 + 1.0/3.0), 1.0e-6);
        assertEquals(ixy, 4.0, 1.0e-6);
    }
}

That cross-product is wrong, which just proves the unit test is incomplete. I notice that you have no polygons with two vertices that don't lie on either axis. — Ben Voigt, Jul 25 '10 at 14:55
Egg on my face - the cross product was indeed incorrect - now corrected. The formulas for Ix, Iy, and Ixy are what you really want, because moment of inertia is a second order tensor. One number isn't enough. — duffymo, Jul 25 '10 at 15:28

score 0 · Answer 3 · edited May 23 '17 at 12:01

0

For reference, here's a mutable 2D org.gcs.kinetic.Vector implementation and a more versatile, immutable org.jscience.mathematics.vector implementation. This article on Calculating a 2D Vector’s Cross Product is helpful, too.

edited May 23 '17 at 12:01

Community

1
1

answered Jul 25 '10 at 15:35

trashgod

203,806
29
246
1,045

score 0 · Answer 4 · answered Nov 07 '10 at 14:44

0

I did it with Tesselation. And take the MOI's all together.

answered Nov 07 '10 at 14:44

Martijn Courteaux

67,591
47
198
287

Calculating the Moment Of Inertia for a concave 2D polygon relative to its orgin

4 Answers4