This is in branch Rambetter-math-fix-experiments.

- Adding new function ChopWindingInPlaceAccu() in polylib.c.  This function is
just like ChopWindingInPlace() except that the winding input has the higher
vect_accu_t resolution.  I actually did a deep massage of this code as well
(went over it with a fine-toothed comb).  I still need to examine it at least
one more time when I'm fresh, and fix an issue related to choosing the value
for maxpts.

- In mathlib.h, defining VEC_SMALLEST_EPSILON and VEC_ACCU_SMALLEST_EPSILON
constants that will help determine suitable epsilons to use for certain
operations.  Comments are provided detailing what these constants do.

- In mathlib.h, defining DotProductAccu() and VectorCopyAccu().

- Small comment in brush.c that has concerns about a certain use of
BaseWindingForPlane() in WriteBSPBrushMap().

State of this code is: compiles on Linux.  May not compile on Windows.  Not
tested (new functionality not used yet in any real code, only functions exist).


git-svn-id: svn://svn.icculus.org/gtkradiant/GtkRadiant/branches/Rambetter-math-fix-experiments@382 8a3a26a2-13c4-0310-b231-cf6edde360e5

libs/mathlib.h

@@ -24,6 +24,7 @@ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 
 // mathlib.h
 #include <math.h>
+#include <float.h>
 
 #include "bytebool.h"
 
@@ -37,6 +38,13 @@ typedef vec_t vec3_t[3];
 typedef vec_t vec5_t[5];
 typedef vec_t vec4_t[4];
 
+// Smallest positive value such that 1.0 + VEC_SMALLEST_EPSILON != 1.0
+// In the case of 32 bits (which is the case), it's 0.00000011921.
+// Don't forget that your epsilons should depend on the possible range of values,
+// because for example 1000.0 + VEC_SMALLEST_EPSILON will almost certainly be
+// equal to 1000.0.
+#define VEC_SMALLEST_EPSILON FLT_EPSILON
+
 #define	SIDE_FRONT		0
 #define	SIDE_ON			2
 #define	SIDE_BACK		1
@@ -308,13 +316,23 @@ vec_t ray_intersect_triangle(const ray_t *ray, qboolean bCullBack, const vec3_t
 typedef double vec_accu_t;
 typedef vec_accu_t vec3_accu_t[3];
 
+// Smallest positive value such that 1.0 + VEC_ACCU_SMALLEST_EPSILON != 1.0
+// In the case of 64 bits (which is the case), it's 0.00000000000000022204.
+// Don't forget that your epsilons should depend on the possible range of values,
+// because for example 1000.0 + VEC_ACCU_SMALLEST_EPSILON will almost certainly
+// be equal to 1000.0.
+#define VEC_ACCU_SMALLEST_EPSILON DBL_EPSILON
+
 // TODO: I have a feeling it may be safer to break these function out into actual functions
 // in order to avoid accidental loss of precision.  For example, say you call
 // VectorScaleAccu(vec3_t, vec_t, vec3_accu_t).  The scale would take place in 32 bit land
 // and the result would be cast to 64 bit, which would cause total loss of precision when
-// scaling by a large factor.
+// scaling by a large factor.  That is why I'm copying and pasting the above #defines here,
+// so that we can break these out into functions if need be.
+#define DotProductAccu(x, y) ((x)[0] * (y)[0] + (x)[1] * (y)[1] + (x)[2] * (y)[2])
 #define VectorSubtractAccu(a, b, c) ((c)[0] = (a)[0] - (b)[0], (c)[1] = (a)[1] - (b)[1], (c)[2] = (a)[2] - (b)[2])
 #define VectorAddAccu(a, b, c) ((c)[0] = (a)[0] + (b)[0], (c)[1] = (a)[1] + (b)[1], (c)[2] = (a)[2] + (b)[2])
+#define VectorCopyAccu(a, b) ((b)[0] = (a)[0], (b)[1] = (a)[1], (b)[2] = (a)[2])
 #define VectorScaleAccu(a, b, c) ((c)[0] = (b) * (a)[0], (c)[1] = (b) * (a)[1], (c)[2] = (b) * (a)[2])
 #define CrossProductAccu(a, b, c) ((c)[0] = (a)[1] * (b)[2] - (a)[2] * (b)[1], (c)[1] = (a)[2] * (b)[0] - (a)[0] * (b)[2], (c)[2] = (a)[0] * (b)[1] - (a)[1] * (b)[0])
 

tools/quake3/common/polylib.c

@@ -592,6 +592,138 @@ void	ClipWindingEpsilon (winding_t *in, vec3_t normal, vec_t dist,
 }
 
 
+/*
+=============
+ChopWindingInPlaceAccu
+=============
+*/
+void ChopWindingInPlaceAccu(winding_accu_t **inout, vec3_t normal, vec_t dist, vec_t crudeEpsilon)
+{
+	vec_accu_t	fineEpsilon;
+	winding_accu_t	*in;
+	int		counts[3];
+	int		i, j;
+	vec_accu_t	dists[MAX_POINTS_ON_WINDING + 1];
+	int		sides[MAX_POINTS_ON_WINDING + 1];
+	int		maxpts;
+	winding_accu_t	*f;
+	vec_accu_t	*p1, *p2;
+	vec_accu_t	w;
+	vec3_accu_t	mid;
+
+	// We require at least a very small epsilon.  It's a good idea for several reasons.
+	// First, we will be dividing by a potentially very small distance below.  We don't
+	// want that distance to be too small; otherwise, things "blow up" with little accuracy
+	// due to the division.  (After a second look, the value w below is in range (0,1), but
+	// graininess problem remains.)  Second, Having minimum epsilon also prevents the following
+	// situation.  Say for example we have a perfect octagon defined by the input winding.
+	// Say our chopping plane (defined by normal and dist) is essentially the same plane
+	// that the octagon is sitting on.  Well, due to rounding errors, it may be that point
+	// 1 of the octagon might be in front, point 2 might be in back, point 3 might be in
+	// front, point 4 might be in back, and so on.  So we could end up with a very ugly-
+	// looking chopped winding, and this might be undesirable, and would at least lead to
+	// a possible exhaustion of MAX_POINTS_ON_WINDING.  It's better to assume that points
+	// very very close to the plane ore on the plane, using an infinitesimal epsilon amount.
+
+	// Now, the original ChopWindingInPlace() function used a vec_t-based winding_t.
+	// So this minimum epsilon is quite similar to casting the higher resolution numbers to
+	// the lower resolution and comparing them in the lower resolution mode.  We explicitly
+	// choose the minimum epsilon as something around the vec_t epsilon of one because we
+	// want the resolution of vec_accu_t to have a large resolution around the epsilon.
+	// Some of that leftover resolution even goes away after we scale to points far away.
+
+	static const vec_accu_t smallestEpsilonAllowed = ((vec_accu_t) VEC_SMALLEST_EPSILON) * 0.5;
+	if (crudeEpsilon < smallestEpsilonAllowed)	fineEpsilon = smallestEpsilonAllowed;
+	else						fineEpsilon = (vec_accu_t) crudeEpsilon;
+
+	in = *inout;
+	counts[0] = counts[1] = counts[2] = 0;
+
+	for (i = 0; i < in->numpoints; i++)
+	{
+		dists[i] = DotProductAccu(in->p[i], normal) - dist;
+		if (dists[i] > fineEpsilon)		sides[i] = SIDE_FRONT;
+		else if (dists[i] < -fineEpsilon)	sides[i] = SIDE_BACK;
+		else					sides[i] = SIDE_ON;
+		counts[sides[i]]++;
+	}
+	sides[i] = sides[0];
+	dists[i] = dists[0];
+	
+	if (!counts[SIDE_FRONT]) {
+		FreeWindingAccu(in);
+		*inout = NULL;
+		return;
+	}
+	if (!counts[SIDE_BACK]) {
+		return; // Winding is unmodified.
+	}
+
+	// TODO: Why the hell are we using this number?  The original comment says something along
+	// the lines of "Can't use counts[0] + 2 because of fp grouping errors", and I'm not
+	// entirely sure what that means.  Is it a problem later on in how the caller of this function
+	// uses the returned winding afterwards?  Is it a problem in that the algorithm below will
+	// bump into the max in some cases?  Anyhow, investigate and fix this number to be a correct
+	// and predictable number.  We can also dynamically grow inside of this function in case we
+	// do hit the max for some strange reason.
+	maxpts = in->numpoints + 4;
+
+	f = AllocWindingAccu(maxpts);
+
+	for (i = 0; i < in->numpoints; i++)
+	{
+		p1 = in->p[i];
+
+		if (sides[i] == SIDE_ON)
+		{
+			if (f->numpoints >= maxpts)
+				Error("ChopWindingInPlaceAccu: points exceeded estimate");
+			if (f->numpoints >= MAX_POINTS_ON_WINDING)
+				Error("ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING");
+			VectorCopyAccu(p1, f->p[f->numpoints]);
+			f->numpoints++;
+			continue;
+		}
+		if (sides[i] == SIDE_FRONT)
+		{
+			if (f->numpoints >= maxpts)
+				Error("ChopWindingInPlaceAccu: points exceeded estimate");
+			if (f->numpoints >= MAX_POINTS_ON_WINDING)
+				Error("ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING");
+			VectorCopyAccu(p1, f->p[f->numpoints]);
+			f->numpoints++;
+		}
+		if (sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i])
+		{
+			continue;
+		}
+			
+		// Generate a split point.
+		p2 = in->p[((i + 1) == in->numpoints) ? 0 : (i + 1)];
+
+		// The divisor's absolute value is greater than the dividend's absolute value.
+		// w is in the range (0,1).
+		w = dists[i] / (dists[i] - dists[i + 1]);
+
+		for (j = 0; j < 3; j++)
+		{
+			// Avoid round-off error when possible.  Check axis-aligned normal.
+			if (normal[j] == 1)		mid[j] = dist;
+			else if (normal[j] == -1)	mid[j] = -dist;
+			else				mid[j] = p1[j] + (w * (p2[j] - p1[j]));
+		}
+		if (f->numpoints >= maxpts)
+			Error("ChopWindingInPlaceAccu: points exceeded estimate");
+		if (f->numpoints >= MAX_POINTS_ON_WINDING)
+			Error("ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING");
+		VectorCopyAccu(mid, f->p[f->numpoints]);
+		f->numpoints++;
+	}
+
+	FreeWindingAccu(in);
+	*inout = f;
+}
+
 /*
 =============
 ChopWindingInPlace

tools/quake3/common/polylib.h

@@ -69,3 +69,4 @@ typedef struct
 } winding_accu_t;
 
 winding_accu_t	*BaseWindingForPlaneAccu(vec3_t normal, vec_t dist);
+void	ChopWindingInPlaceAccu(winding_accu_t **w, vec3_t normal, vec_t dist, vec_t epsilon);

tools/quake3/q3map2/brush.c

@@ -506,6 +506,8 @@ void WriteBSPBrushMap( char *name, brush_t *list )
 		fprintf (f, "{\n");
 		for (i=0,s=list->sides ; i<list->numsides ; i++,s++)
 		{
+			// TODO: See if we can use a smaller winding to prevent resolution loss.
+			// Is WriteBSPBrushMap() used only to decompile maps?
 			w = BaseWindingForPlane (mapplanes[s->planenum].normal, mapplanes[s->planenum].dist);
 
 			fprintf (f,"( %i %i %i ) ", (int)w->p[0][0], (int)w->p[0][1], (int)w->p[0][2]);