Get the Surface Area of a Polyhedron (3D object)

问题:

I have a 3D surface, (think about the xy plane). The plane can be slanted. (think about a slope road).

Given a list of 3D coordinates that define the surface(Point3D1X, Point3D1Y, Point3D1Z, Point3D12X, Point3D2Y, Point3D2Z, Point3D3X, Point3D3Y, Point3D3Z, and so on), how to calculate the area of the surface?

Note that my question here is analogous to finding area in 2D plane. In 2D plane we have a list of points that defines a polygon, and using this list of points we can find the area of the polygon. Now assuming that all these points have z values in such a way that they are elevated in 3D to form a surface. My question is how to find the area of that 3D surface?

回答1:

Since you say it's a polyhedron, stacker's link (http://softsurfer.com/Archive/algorithm_0101/algorithm_0101.htm) is applicable.

Here's my approximate C# translation of the C code for your situation:

// NOTE: The original code contained the following notice:
// ---------------------------------------
// Copyright 2000 softSurfer, 2012 Dan Sunday
// This code may be freely used and modified for any purpose
// providing that this copyright notice is included with it.
// iSurfer.org makes no warranty for this code, and cannot be held
// liable for any real or imagined damage resulting from its use.
// Users of this code must verify correctness for their application.
// ---------------------------------------
// area3D_Polygon(): computes the area of a 3D planar polygon
//    Input:  int n = the number of vertices in the polygon
//            Point[] V = an array of n+2 vertices in a plane
//                       with V[n]=V[0] and V[n+1]=V[1]
//            Point N = unit normal vector of the polygon's plane
//    Return: the (float) area of the polygon
static float
area3D_Polygon( int n, Point3D[] V, Point3D N )
{
    float area = 0;
    float an, ax, ay, az;  // abs value of normal and its coords
    int   coord;           // coord to ignore: 1=x, 2=y, 3=z
    int   i, j, k;         // loop indices

    // select largest abs coordinate to ignore for projection
    ax = (N.x>0 ? N.x : -N.x);     // abs x-coord
    ay = (N.y>0 ? N.y : -N.y);     // abs y-coord
    az = (N.z>0 ? N.z : -N.z);     // abs z-coord

    coord = 3;                     // ignore z-coord
    if (ax > ay) {
        if (ax > az) coord = 1;    // ignore x-coord
    }
    else if (ay > az) coord = 2;   // ignore y-coord

    // compute area of the 2D projection
    for (i=1, j=2, k=0; i<=n; i++, j++, k++)
        switch (coord) {
        case 1:
            area += (V[i].y * (V[j].z - V[k].z));
            continue;
        case 2:
            area += (V[i].x * (V[j].z - V[k].z));
            continue;
        case 3:
            area += (V[i].x * (V[j].y - V[k].y));
            continue;
        }

    // scale to get area before projection
    an = Math.Sqrt( ax*ax + ay*ay + az*az);  // length of normal vector
    switch (coord) {
    case 1:
        area *= (an / (2*ax));
        break;
    case 2:
        area *= (an / (2*ay));
        break;
    case 3:
        area *= (an / (2*az));
        break;
    }
    return area;
}

回答2:

I upvoted a few answers which I think are correct. But I think the simplest way to do it-- regardless of whether it's in 2D or 3D, is to use the following formula:

area = sum(V(i+1)XV(i))/2;

Where X is the vector cross.

The code to do this is:

    public double Area(List<Point3D> PtList)
    {

        int nPts = PtList.Count;
        Point3D a;
        int j = 0;

        for (int i = 0; i < nPts; ++i)
        {
            j = (i + 1) % nPts;
            a += Point3D.Cross(PtList[i], PtList[j]);
        }
        a /= 2;
        return Point3D.Distance(a,default(Point3D));
    }

    public static Point3D Cross(Point3D v0, Point3D v1)
    {
        return new Point3D(v0.Y * v1.Z - v0.Z * v1.Y,
            v0.Z * v1.X - v0.X * v1.Z,
            v0.X * v1.Y - v0.Y * v1.X);
    }

Note that the solution doesn't depend on projection to x-plane, which I think is clunky.

What do you think?

回答3:

Do you mean Area of 3D planer polygons?

1. http://softsurfer.com/Archive/algorithm_0101/algorithm_0101.htm
2. http://local.wasp.uwa.edu.au/~pbourke/geometry/area3d/
3. http://thebuildingcoder.typepad.com/blog/2008/12/3d-polygon-areas.html

回答4:

I don't know about optimising this method (I've not done it in code before), but the way to approach it mathematically is to split your shape into triangles, the area of which is then easily calculated and summed. (Remember: area of a triangle is width * height * 0.5 - you'll need to calculate the height of non-right-angled triangles.)

Doing these things in 3D generally means one more calculation is needed at each stage. For example, in 2D, the distance between 2 points (the length of a side of your shape) is calculated something like this (pseduocode because I don't have VS on this machine):

double DistanceBetween(Point a, Point b)
{
   double dx = a.x - b.x;
   double dy = a.y - b.y;
   return SquareRoot(dx*dx + dy*dy);
}

In three dimensions that becomes:

double DistanceBetween(Point3d a, Point3d b)
{
   double dx = a.x - b.x;
   double dy = a.y - b.y;
   double dz = a.z - b.z;
   return SquareRoot(dx*dx + dy*dy + dz*dz);
}

Splitting a shape into arbitrary triangles just involves picking any three adjacent vertices at a time, until your're down to your last three.

回答5:

You can derive the solution in terms of the 2D solution.

Consider the polygon mad up from a heap of smaller triangles.

Project each triangle back to the XY plane. You can show that the area of the orignal triangle is 1/(n.k) times the area of the projected triangle. (Here n is the unit normal to the plane containing the polygon, and k is the unit vector in the z direction)

So the total area of the original is 1/(n.k) times the area of the polygon projected into the XY plane. Which you can work out using your existing 2D formula.

You can calculate n by (e1 x e2 ) / || e1 x e2 || where e1 and e2 are any 2 non-parallel edges of your polygon.

Of course you may get better (more accurate) results by projecting into the XZ or YZ plane.. you should pick the one with the normal closest to that of your plane.

回答6:

Another solution that won't require you to create a mesh of polygons is to do a contour integral around the perimeter. You use Green's theorem to convert the area integral into a contour integral, then use something simple like Gauss quadrature to integrate and sum each contribution. You do have to have a definition of the perimeter.

This process can work with 2D shapes that have holes as well. You just have to define a cut that runs from the outer perimeter to the hole, integrate around the hole, and then go back out to the perimeter.

回答7:

@Graviton I can't comment on the answer above so I will submit a new one.

This might be my unfamiliarity with c# syntax, but I believe your answer is missing the dot product with the unit normal vector. The formula should be:

area = n.sum( V(i+1) x V(i) )/2;

where n refers to unit normal vector to the plane, . to dot product and x cross product.

The normal can be calculated using any 3 vectors of the polygon:

n = (V1-V0)x(V2-V0)/magnitude((V1-V0)x(V2-V0))

Here's a javascript implementation using the Vector.js lib:

  function getArea (vecs) {
    var area = 0;
    var vecs = [];
    var j = 0;
    var a = new Vector(0,0,0);

    for (var i = 0; i < vecs.length; i++) {
      j = (i + 1) % vecs.length;
      a = a.add( vecs[i].cross(vecs[j]) );
    }
    a = a.divide(2);
    var v1 = vecs[1].subtract(vecs[0]);
    var v2 = vecs[2].subtract(vecs[0]);
    var normal = v1.cross(v2);
    normal = normal.unit();
    // area = a.length()/10000; // convert to m2
    area = (normal.dot(a))/10000;
    return area;
  };