views:

400

answers:

6

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( Point3D1, Point3D2, Point3D3, 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 plane/surface. My question is how to find the area of that 3D plane/surface?

+1  A: 

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.

Michael Anderson
How are you going to generalize the formula to a list of points, inside of just 3 points?
Ngu Soon Hui
+1  A: 

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
+7  A: 

Since you say it's a polygon, 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:

// 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;
}
Gabe
+1 For the code:)
TheMachineCharmer
@The Machine Charmer you forgot to upvote ;-) did it
stacker
@stacker Done now ;)
TheMachineCharmer
+1  A: 

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?

Ngu Soon Hui
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
@gabe: The solution you proposed requires projection, which I don't like. Otherwise it's a good one.
Ngu Soon Hui
+1  A: 

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