计算平面上闭合多边形的面积



我正在尝试计算位于平面上的多边形的面积(形成不相交闭合形状的共面点的集合),我知道一种方法可以计算二维(而不是三维)不规则(或任何)多边形的面积。我的解决方案是旋转平面,使其在z方向上的法线为0(这样我就可以像对待2D一样对待它),然后运行2D面积函数。

问题是,我不知道如何实际确定旋转轴和使平面在Z轴上变平的量。我通过我能找到的最简单的三维旋转方法进行旋转:旋转矩阵。那么,考虑到我试图使用旋转矩阵来进行旋转,我如何计算出旋转平面的角度,使其与另一个向量的方向相同?实际上,我不太懂微积分或欧几里得几何,所以无论哪种解决方案需要我自学最少,都是理想的解决方案。有更好的方法吗?

下面是我的尝试,它甚至还没有达到平面在Z轴上。这是我的"Surface"类的一个实例方法,它是我"Plane"类的衍生物,并且有一个形成闭合多边形的共面点(IntersectPoints)阵列。

public virtual double GetArea()
{
Vector zUnit = new Vector(0, 0, 1); //vector perprendicualr to z
Vector nUnit = _normal.AsUnitVector();
Surface tempSurface = null;
double result = 0;
if (nUnit != zUnit && zUnit.Dot(nUnit) != 0) //0 = perprendicular to z
{
tempSurface = (Surface)Clone();
double xAxisAngle = Vector.GetAxisAngle(nUnit, zUnit, Physics.Formulae.Axes.X);
double yAxisAngle = Vector.GetAxisAngle(nUnit, zUnit, Physics.Formulae.Axes.Y);
double rotationAngle = Vector.GetAxisAngle(nUnit, zUnit, Physics.Formulae.Axes.Z);
tempSurface.Rotate(xAxisAngle, yAxisAngle, rotationAngle); //rotating plane so that it is flat on the Z axis
}
else
{
tempSurface = this;
}
for (int x = 0; x < tempSurface.IntersectPoints.Count; x++) //doing a cross sum of each point
{
Point curPoint = tempSurface.IntersectPoints[x];
Point nextPoint;
if (x == tempSurface.IntersectPoints.Count - 1)
{
nextPoint = tempSurface.IntersectPoints[0];
}
else
{
nextPoint = tempSurface.IntersectPoints[x + 1];
}
double cross1 = curPoint.X * nextPoint.Y;
double cross2 = curPoint.Y * nextPoint.X;
result += (cross1 - cross2); //add the cross sum of each set of points to the result
}
return Math.Abs(result / 2); //divide cross sum by 2 and take its absolute value to get the area.
}

下面是我的核心旋转和获取轴角度的方法:

private Vector Rotate(double degrees, int axis)
{
if (degrees <= 0) return this;
if (axis < 0 || axis > 2) return this;
degrees = degrees * (Math.PI / 180); //convert to radians
double sin = Math.Sin(degrees);
double cos = Math.Cos(degrees);
double[][] matrix = new double[3][]; 
//normalizing really small numbers to actually be zero
if (Math.Abs(sin) < 0.00000001) 
{
sin = 0;
}
if (Math.Abs(cos) < 0.0000001)
{
cos = 0;
}
//getting our rotation matrix
switch (axis)
{
case 0: //x axis
matrix = new double[][] 
{ 
new double[] {1, 0, 0},
new double[] {0, cos, sin * -1},
new double[] {0, sin, cos}
};
break;
case 1: //y axis
matrix = new double[][] 
{ 
new double[] {cos, 0, sin},
new double[] {0, 1, 0},
new double[] {sin * -1, 0, cos}
};
break;
case 2: //z axis
matrix = new double[][] 
{ 
new double[] {cos, sin * -1, 0},
new double[] {sin, cos, 0},
new double[] {0, 0, 1}
};
break;
default:
return this;
}
return Physics.Formulae.Matrix.MatrixByVector(this, matrix);
}
public static double GetAxisAngle(Point a, Point b, Axes axis, bool inDegrees = true)
{ //pretty sure this doesnt actually work
double distance = GetDistance(a, b);
double difference;
switch (axis)
{
case Axes.X:
difference = b.X - a.X;
break;
case Axes.Y:
difference = b.Y - a.Y;
break;
case Axes.Z :
difference = b.Z - a.Z;
break;
default:
difference = 0;
break;
}
double result = Math.Acos(difference / distance);
if (inDegrees == true)
{
return result * 57.2957; //57.2957 degrees = 1 radian
}
else
{
return result;
}
}

实现这一点的一种稳健方法是对每条边的顶点的叉积求和。如果顶点是共面的,这将生成平面的法线,其长度是闭合多边形面积的2倍。

请注意,此方法与问题中链接的二维方法非常相似,后者实际上计算三维叉积的二维等效值,对所有边求和,然后除以2。

Vector normal = points[count-1].cross(points[0]);
for(int i=1; i<count; ++i) {
normal += points[i-1].cross(points[i]);
}
double area = normal.length() * 0.5;

这种方法的优点:

  • 如果你的顶点只是近似平面的,它仍然给出了正确的答案
  • 它不取决于平面的角度
  • 事实上,你根本不需要处理这个角度
  • 如果希望知道平面方向,则已经获得法线

一个可能的困难是:如果你的多边形很小,离原点很远,你可能会遇到浮点精度问题。如果可能出现这种情况,您应该首先平移所有顶点,使其中一个位于原点,如下所示:

Vector normal(0,0,0);
Vector origin = points[count-1]; 
for(int i=1; i<count-1; ++i) {
normal += (points[i-1]-origin).cross(points[i]-origin);
}
double area = normal.length() * 0.5;

您不需要旋转平面(或所有点)。例如,只需使用GetArea函数计算多边形投影到Z平面的面积(如果它不垂直于多边形平面),并将结果除以多边形平面的余弦-Z平面角度-它等于zUnit和nUnit的标量乘积(我建议nUnit是多边形平面的法向量)

TrueArea = GetArea() / zUnit.Dot(nUnit)

最新更新