Calculating the Area of a Closed Polygon on a Plane

Question

I'm attempting to calculate the area of a polygon that lies on a plane (a collection co-planar points forming a non-intersecting closed shape), and I know a method that can calculate the area of an irregular (or any) polygon in two dimensions - but not three. My solution is to rotate the plane so that it's normal is 0 in the z direction (so I can treat it like it's 2D) and then run the 2D area function.

The problem is I have NO idea how to actually determine the rotation axes and amounts to flatten a plane on it's Z-axis. I do my rotation through the easiest method I could find for 3 dimensional rotation: Rotation Matrices. So, given that I'm trying to use rotation matrices to do my rotation, how do I figure out the angles to rotate my plane by to be oriented in the same direction as another vector? I don't actually know much calculus or Euclidean geometry, so whichever solution requires me to teach myself the least of both is the ideal solution. Is there a better way?

Here's my attempt below, which doesn't even come close to getting the plane flat on the Z axis. This is an instance method of my "Surface" class, which is a derivative of my "Plane" class, and has an array of co-planar points (IntersectPoints) forming a closed polygon.

 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.
    }

And here are my core rotation and get axis angle methods:

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

Answer 1

A robust way to do this is to do a sum of the cross-products of the vertices of each edge. If your vertices are co-planar, this will produce a normal to the plane, whose length is 2 times the area of the closed polygon.

Note that this method is very similar to the 2D method linked in your question, which actually calculates a 2D equivalent of the 3D cross-product, summed for all edges, then divides by 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;

Advantages of this method:

• If your vertices are only approximately planar, it still gives the right answer
• It doesn't depend on the angle of the plane.
• In fact you don't need to deal with the angle at all.
• If you want to know the plane orientation, you've got the normal already.

One possible difficulty: if your polygon is very small, and a long way away from the origin, you can get floating point precision problems. If that case is likely to arise, you should first translate all of your vertices so that one is at the origin, like so:

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;

Answer 2

You need not to rotate the plane (or all points). Just calculate an area of polygon projection to Z-plane (if it is not perpendicular to polygon plane), for example, with you GetArea function, and divide result by cosinus of Poly-plane - Z-plane angle - it is equal to scalar product of zUnit and nUnit (I suggest that nUnit is normal vector to polygon plane)

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