12

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?

Graviton
  • 81,782
  • 146
  • 424
  • 602

7 Answers7

11

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;
}
Peter O.
  • 32,158
  • 14
  • 82
  • 96
Gabe
  • 84,912
  • 12
  • 139
  • 238
4

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) × V(i))/2;

Where × 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?

StayOnTarget
  • 11,743
  • 10
  • 52
  • 81
Graviton
  • 81,782
  • 146
  • 424
  • 602
  • Other solutions (like I posted) do fewer math operations (1 multiply, 1 subtract) per point. Yours does 6 multiplies and 3 subtracts per point. If you don't have too many points or don't have to do it often enough, use the method that seems simpler. Otherwise if one is fast enough to notice the difference, use the faster one. – Gabe Feb 28 '10 at 11:03
  • @gabe: The solution you proposed requires projection, which I don't like. Otherwise it's a good one. – Graviton Feb 28 '10 at 11:40
  • @Graviton: Are you certain that this is correct? I believe it does not hold true for the general case. You are summing the area of each triangle along a list of points. Are you assuming that this list of points actually constitutes a convex polyon? Either way I do not see how this can result in the correct area of a surface as (1) you will have overlapping triangles and (2) you will miss areas. – Chris Dec 07 '12 at 10:27
  • If you are assuming an ordered list of points resulting in a convex polygon, I assume you would want to triangulate correctly - in this case, one should not use PtList[i], PtList[j], but keep one (e.g. the first) point constant. Meaning: you don't need to use i and calculating the local area would look like this: Point3D.Cross(PtList[0], PtList[j]); – Chris Dec 07 '12 at 10:36
  • @Chris, this is a direct extension of the 2D Area formula. And in 2D, the equivalent formula doesnt' really care whether you are a convex or concave. – Graviton Dec 07 '12 at 11:12
1

@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;
  };
winkerVSbecks
  • 1,173
  • 1
  • 10
  • 24
1

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

Consider the polygon made 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.

StayOnTarget
  • 11,743
  • 10
  • 52
  • 81
Michael Anderson
  • 70,661
  • 7
  • 134
  • 187
1

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.

Dan Puzey
  • 33,626
  • 4
  • 73
  • 96
1

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.

duffymo
  • 305,152
  • 44
  • 369
  • 561