1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
#pragma once

// Functions from paper entitled "Fast Polygon Area and Newell Normal Computation",
// Daniel Sunday, journal of graphics tools, 7(2):9-13, 2002.
// http://www.acm.org/jgt/papers/Sunday02/FastArea.html

#define PREV(i, N) ((i + N-1) % N)
#define NEXT(i, N) ((i + 1) % N)

// Input: 2D polygon
static inline double findArea(int n, double *x, double *y)        
{
	// Assume that the 2D polygon's vertex coordinates are
	// stored in (x[0], y[0]), ..., (x[n-1], y[n-1]),
	// with no assumed vertex replication at the end.

	// Initialize sum with the boundary vertex terms
	double sum = x[0] * (y[1] - y[n-1]) + x[n-1] * (y[0] - y[n-2]);

	for (int i=1; i < n-1; i++) {
		sum += x[i] * ( y[i+1] - y[i-1] );
	}

	return (sum / 2.0);
}

// return the signed area of a 3D planar polygon (given normal vector)
// 3D planar polygon, and plane normal
static inline double findArea3D(int n, double *x, double *y, double *z, double nx, double ny, double nz) 
{
	// select largest normal coordinate to ignore for projection
	double ax = (nx>0 ? nx : -nx);	// abs nx
	double ay = (ny>0 ? ny : -ny);	// abs ny
	double az = (nz>0 ? nz : -nz);	// abs nz
	double len = sqrt(nx*nx + ny*ny + nz*nz); // length of normal

	if (ax > ay) {
		if (ax > az)			       // ignore x-coord
			return findArea(n, y, z) * (len / nx);
	}
	else if (ay > az)			       // ignore y-coord
		return findArea(n, z, x) * (len / ny);
	return findArea(n, x, y) * (len / nz); // ignore z-coord
}

// output the approximate unit normal of a 3D nearly planar polygon
// return the area of the polygon
// Inputs: 3D polygon, output unit normal
static inline double findNormal3D(int n, double *x, double *y, double *z, double *nx, double *ny, double *nz) 
{
	// get the Newell normal
	double nwx = findArea(n, y, z);
	double nwy = findArea(n, z, x);
	double nwz = findArea(n, x, y);
	// get length of the Newell normal
	double nlen = sqrt( nwx*nwx + nwy*nwy + nwz*nwz );
	// compute the unit normal
	*nx = nwx / nlen;
	*ny = nwy / nlen;
	*nz = nwz / nlen;
	return nlen;    // area of polygon = length of Newell normal
}

static inline Vector3 centerOfPoints(const std::vector< Vector3 > & points)
{
	Vector3 center;

	for(std::vector<Vector3>::const_iterator it = points.begin(); it != points.end(); it++)
		center += *it;

	center /= points.size();

	return center;
}

static inline double findArea3D(const std::vector<Vector3>& points, const Vector3& n)
{
	int N = points.size();

	std::vector<double> x(N), y(N), z(N);

	for(int i = 0; i < N; i++)
	{
		x[i] = points[i][0];
		y[i] = points[i][1];
		z[i] = points[i][2];
	}

	double * X = &x.front();
	double * Y = &y.front();
	double * Z = &z.front();

	double area = findArea3D(N, X, Y, Z, n[0], n[1], n[2]);

	return area;
}

static inline double findNormal3D(const std::vector<Vector3>& points, Vector3& n)
{
	int N = points.size();

	std::vector<double> x(N), y(N), z(N);

	for(int i = 0; i < N; i++)
	{
		x[i] = points[i][0];
		y[i] = points[i][1];
		z[i] = points[i][2];
	}

	double * X = &x.front();
	double * Y = &y.front();
	double * Z = &z.front();

	return findNormal3D(N, X, Y, Z, &n[0], &n[1], &n[2]);
}

static inline double signedArea(const std::vector<Vector3>& points, const Vector3& n, const Vector3& center)
{
	int N = points.size();
	Vector3 sumVec;

	for(int i = 0; i < N; i++)
	{
		int j = NEXT(i, N);
		sumVec += cross(Vector3(points[i]-center),Vector3(points[j]-center));
	}

	return dot(sumVec, n);
}