MeshKit  1.0
FastArea.hpp
Go to the documentation of this file.
```00001 /*
00002  * FastArea.c++
00003  *
00004  * From the paper:
00005  *
00006  *      Daniel Sunday
00007  *      "Fast Polygon Area and Newell Normal Computation"
00008  *      journal of graphics tools, 7(2):9-13, 2002
00009  *
00010  */
00011
00012 // assume vertex coordinates are in arrays x[], y[], and z[]
00013
00014 // with room to duplicate the first two vertices at the end
00015
00016
00017 // return the signed area of a 2D polygon
00018
00019 #include <math.h>
00020 #include <iostream>
00021
00022 using namespace std;
00023
00024 inline double
00025 PolygonArea(int n, double *x, double *y)             // 2D polygon
00026
00027 {
00028     // guarantee the first two vertices are also at array end
00029
00030     x[n] = x[0];    y[n] = y[0];
00031     x[n+1] = x[1];  y[n+1] = y[1];
00032
00033     double sum = 0.0;
00034     double *xptr = x+1, *ylow = y, *yhigh = y+2;
00035     for (int i=1; i <= n; i++) {
00036         sum += (*xptr++) * ( (*yhigh++) - (*ylow++) );
00037     }
00038
00039     return 0.5*sum ;
00040 }
00041
00042 // return the signed area of a 3D planar polygon (given normal vector)
00043
00044 inline double
00045 PolygonArea3D(int n, double *x, double *y, double *z,  // 3D planar polygon
00046
00047            double nx, double ny, double nz)         // and plane normal
00048
00049 {
00050     // select largest normal coordinate to ignore for projection
00051     // CSV: Original Code seeems to have bug: Need to verify it.
00052
00053     double ax = (nx>0 ? nx : -nx);      // abs nx
00054
00055     double ay = (ny>0 ? ny : -ny);      // abs ny
00056
00057     double az = (nz>0 ? nz : -nz);      // abs nz
00058
00059     double len = sqrt(nx*nx + ny*ny + nz*nz); // length of normal
00060
00061
00062     if (ax > ay) {
00063         if (ax > az)                           // ignore x-coord
00064
00065 // CSV   return PolygonArea(n, y, z) * (len / nx);
00066            return PolygonArea(n, y, z) * (len / ax);
00067     }
00068     else if (ay > az)                          // ignore y-coord
00069
00070 // CSV   return PolygonArea(n, z, x) * (len / ny);
00071      return PolygonArea(n, z, x) * (len / ay);
00072
00073 //CSV   return PolygonArea(n, x, y) * (len / nz); // ignore z-coord
00074    return PolygonArea(n, x, y) * (len / az); // ignore z-coord
00075
00076 }
00077
00078 // output the approximate unit normal of a 3D nearly planar polygon
00079 // return the area of the polygon
00080
00081 inline double
00082 PolygonNormal3D(int n, double *x, double *y, double *z, // 3D polygon
00083              double *nx, double *ny, double *nz) // output unit normal
00084
00085 {
00086     // get the Newell normal
00087
00088     double nwx = PolygonArea(n, y, z);
00089     double nwy = PolygonArea(n, z, x);
00090     double nwz = PolygonArea(n, x, y);
00091
00092     // get length of the Newell normal
00093
00094     double nlen = sqrt( nwx*nwx + nwy*nwy + nwz*nwz );
00095     // compute the unit normal
00096
00097     double multby = 1.0/( double)nlen;
00098
00099     *nx = multby*nwx;
00100     *ny = multby*nwy;
00101     *nz = multby*nwz;
00102
00103     return nlen;    // area of polygon = length of Newell normal
00104 }
00105
00106 inline
00107 double PolygonArea3D( int n, double *x, double *y, double *z)
00108 {
00109     double nx, ny, nz;
00110     PolygonNormal3D( n, x, y, z, &nx, &ny, &nz );
00111
00112     return PolygonArea3D( n, x, y, z, nx, ny, nz );
00113 }
00114
```