计算凹多边形相对于其原点的转动惯量

Question

我希望计算一个(2D)凹多边形的惯性矩。我在网上找到了这个公式。但我不太确定如何解释这个公式... 公式 http://img101.imageshack.us/img101/8141/92175941c14cadeeb956d8f.gif 1) 这个公式正确吗？ 2) 如果是，我的C++转换是否正确？

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

- Martijn Courteaux

---

这实际上不需要对LaTeX的了解。 - Larry Wang

再看一遍，你所说的关于叉积的部分对我来说并没有什么意义。叉积 只适用于三维向量，而你说你正在使用二维向量... - Larry Wang

@Larry Wang：我还没有对上述公式进行过多（或任何）思考，但因为R ^ 2同构于R ^ 3中的平面z = 0，所以您实际上也可以使用叉积来处理二维向量，隐式地将它们的z坐标设置为零。 - Andreas Rejbrand

@Andreas：那是解决那个问题的方案吗？我在转换Larry的代码时遇到了那个问题.... - Martijn Courteaux

你代码中的一些注释有误。叉积在二维和三维空间中总是创建一个向量。只是在二维空间中，它恰好是垂直于二维平面的向量。它绝不是标量。不存在标量和向量的叉积。您可以将向量乘以标量以获得另一个向量，但这不是叉积。 - duffymo

Answer 1

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

编辑：进行了许多小的数学更改

- Larry Wang

---

这实际上不是叉积的长度（大小或二范数），而是叉积的z坐标。任何向量的长度都是非负的。不确定哪个是正确的，但公式似乎要求长度。 - Ben Voigt

而且，函数length的命名是错误的（它是长度的平方）。 - Ben Voigt

Answer 2

我认为你需要做的不仅仅是将公式翻译成代码。你需要确切地理解这个公式的含义。

当你有一个二维多边形时，你可以相对于给定的坐标系计算三个惯性矩：绕x轴的惯性矩、绕y轴的惯性矩和极惯性矩。有一个平行轴定理可以让你从一个坐标系转换到另一个坐标系。

你是否确切知道这个公式适用于哪个惯性矩和坐标系？

以下是一些可能会帮助你的代码，以及一个JUnit测试来证明它的有效性：

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 https://dev59.com/RE_Sa4cB1Zd3GeqP_VDH
     */
    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;
    }
}

这是与之配套的JUnit测试：

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

- duffymo

---

那个叉积是错误的，这证明了单元测试不完整。我注意到你没有多边形有两个顶点不在任何轴上。 - Ben Voigt

此外，你所有的测试用例都是凸多边形。 - Ben Voigt

我很尴尬 - 叉积确实是错误的 - 现在已经更正了。Ix，Iy和Ixy的公式才是你真正想要的，因为惯性矩是一个二阶张量。一个数字是不够的。 - duffymo

Answer 3

参考以下可变的二维org.gcs.kinetic.Vector实现和更灵活、不可变的org.jscience.mathematics.vector实现。此外，计算二维向量叉积的文章也很有帮助。

- trashgod

Answer 4

我用镶嵌完成了它。并将所有的MOI放在一起。

- Martijn Courteaux