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

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

#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

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

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.

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