计算凹二维多边形相对于其原点的惯性矩答案

【问题标题】：Calculating the Moment Of Inertia for a concave 2D polygon relative to its orgin计算凹二维多边形相对于其原点的惯性矩
【发布时间】：2011-03-20 18:04:21
【问题描述】：

我想计算（2D）凹多边形的转动惯量。我在互联网上找到了这个。但我不太清楚如何解释公式...

Formula 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;

马丁

【问题讨论】：

这实际上不需要任何 LaTeX 知识
再看一遍，你说的关于交叉产品的部分对我来说真的没有意义。 Cross product 仅针对 3-D 向量定义，而您说您使用的是 2-D 向量...
@Larry Wang：我对上面的公式没有太多（任何）想法，但是因为 R^2 与 R^3 中的平面 z = 0 同构，所以您实际上可以使用叉积对于 2D 向量也是如此，隐式地将它们的 z 坐标设置为零。
@Andreas：这是解决这个问题的方法吗？我在转换拉里的代码时遇到了这个问题....
@Martijn：正如我所说，我没有仔细研究过这个公式，但正如 Larry Wang 指出的那样，叉积只能在三个维度上计算。但通常的做法是在平面中使用它，将平面假装为三维空间的子空间 z = 0，即在所有平面向量中添加一个等于 0 的第三个 z 坐标。

标签： java c++ math vector 2d

【解决方案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);

编辑：许多小的数学变化

【讨论】：

这实际上不是叉积的长度（大小或二范数），而是叉积的 z 坐标。任何向量的长度都是非负的。不确定哪个是正确的，但公式似乎要求长度。
而且，函数length 命名错误（它是长度的平方）。

【解决方案2】：

我认为您还有更多工作要做，只是将公式转换为代码。您需要准确理解这个公式的含义。

当您有一个 2D 多边形时，您可以计算相对于给定坐标系的三个惯性矩：关于 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 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;
    }
}

这是伴随它的 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);
    }
}

【讨论】：

那个叉积是错误的，只是证明单元测试不完整。我注意到您没有两个顶点不在任一轴上的多边形。
另外，你所有的测试用例都是凸多边形。
我脸上的鸡蛋 - 叉积确实不正确 - 现在已更正。 Ix、Iy 和 Ixy 的公式是您真正想要的，因为转动惯量是二阶张量。一个数字是不够的。

【解决方案3】：

作为参考，这里有一个可变的 2D org.gcs.kinetic.Vector 实现和一个更通用、不可变的 org.jscience.mathematics.vector 实现。 Calculating a 2D Vector’s Cross Product 上的这篇文章也很有帮助。

【讨论】：

【解决方案4】：

我用 Tesselation 做到了。并把 MOI 放在一起。

【讨论】：