// mathlib.h
#include <math.h>
+#include <float.h>
#include "bytebool.h"
typedef vec_t vec5_t[5];
typedef vec_t vec4_t[4];
+// Smallest positive value for vec_t such that 1.0 + VEC_SMALLEST_EPSILON_AROUND_ONE != 1.0.
+// In the case of 32 bit floats (which is almost certainly the case), it's 0.00000011921.
+// Don't forget that your epsilons should depend on the possible range of values,
+// because for example adding VEC_SMALLEST_EPSILON_AROUND_ONE to 1024.0 will have no effect.
+#define VEC_SMALLEST_EPSILON_AROUND_ONE FLT_EPSILON
+
#define SIDE_FRONT 0
#define SIDE_ON 2
#define SIDE_BACK 1
#define Q_rint(in) ((vec_t)floor(in+0.5))
+qboolean VectorIsOnAxis(vec3_t v);
+qboolean VectorIsOnAxialPlane(vec3_t v);
+
vec_t VectorLength(vec3_t v);
void VectorMA( const vec3_t va, vec_t scale, const vec3_t vb, vec3_t vc );
/*! return true if triangle intersects ray... dist = dist from intersection point to ray-origin */
vec_t ray_intersect_triangle(const ray_t *ray, qboolean bCullBack, const vec3_t vert0, const vec3_t vert1, const vec3_t vert2);
+
+////////////////////////////////////////////////////////////////////////////////
+// Below is double-precision math stuff. This was initially needed by the new
+// "base winding" code in q3map2 brush processing in order to fix the famous
+// "disappearing triangles" issue. These definitions can be used wherever extra
+// precision is needed.
+////////////////////////////////////////////////////////////////////////////////
+
+typedef double vec_accu_t;
+typedef vec_accu_t vec3_accu_t[3];
+
+// Smallest positive value for vec_accu_t such that 1.0 + VEC_ACCU_SMALLEST_EPSILON_AROUND_ONE != 1.0.
+// In the case of 64 bit doubles (which is almost certainly the case), it's 0.00000000000000022204.
+// Don't forget that your epsilons should depend on the possible range of values,
+// because for example adding VEC_ACCU_SMALLEST_EPSILON_AROUND_ONE to 1024.0 will have no effect.
+#define VEC_ACCU_SMALLEST_EPSILON_AROUND_ONE DBL_EPSILON
+
+vec_accu_t VectorLengthAccu(const vec3_accu_t v);
+
+// I have a feeling it may be safer to break these #define functions 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.
+//#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])
+//#define Q_rintAccu(in) ((vec_accu_t) floor(in + 0.5))
+
+vec_accu_t DotProductAccu(const vec3_accu_t a, const vec3_accu_t b);
+void VectorSubtractAccu(const vec3_accu_t a, const vec3_accu_t b, vec3_accu_t out);
+void VectorAddAccu(const vec3_accu_t a, const vec3_accu_t b, vec3_accu_t out);
+void VectorCopyAccu(const vec3_accu_t in, vec3_accu_t out);
+void VectorScaleAccu(const vec3_accu_t in, vec_accu_t scaleFactor, vec3_accu_t out);
+void CrossProductAccu(const vec3_accu_t a, const vec3_accu_t b, vec3_accu_t out);
+vec_accu_t Q_rintAccu(vec_accu_t val);
+
+void VectorCopyAccuToRegular(const vec3_accu_t in, vec3_t out);
+void VectorCopyRegularToAccu(const vec3_t in, vec3_accu_t out);
+vec_accu_t VectorNormalizeAccu(const vec3_accu_t in, vec3_accu_t out);
+
#ifdef __cplusplus
}
#endif
vec3_t vec3_origin = {0.0f,0.0f,0.0f};
+/*
+================
+VectorIsOnAxis
+================
+*/
+qboolean VectorIsOnAxis(vec3_t v)
+{
+ int i, zeroComponentCount;
+
+ zeroComponentCount = 0;
+ for (i = 0; i < 3; i++)
+ {
+ if (v[i] == 0.0)
+ {
+ zeroComponentCount++;
+ }
+ }
+
+ if (zeroComponentCount > 1)
+ {
+ // The zero vector will be on axis.
+ return qtrue;
+ }
+
+ return qfalse;
+}
+
+/*
+================
+VectorIsOnAxialPlane
+================
+*/
+qboolean VectorIsOnAxialPlane(vec3_t v)
+{
+ int i;
+
+ for (i = 0; i < 3; i++)
+ {
+ if (v[i] == 0.0)
+ {
+ // The zero vector will be on axial plane.
+ return qtrue;
+ }
+ }
+
+ return qfalse;
+}
+
/*
================
MakeNormalVectors
}
vec_t VectorNormalize( const vec3_t in, vec3_t out ) {
- vec_t length, ilength;
- length = (vec_t)sqrt (in[0]*in[0] + in[1]*in[1] + in[2]*in[2]);
+ // The sqrt() function takes double as an input and returns double as an
+ // output according the the man pages on Debian and on FreeBSD. Therefore,
+ // I don't see a reason why using a double outright (instead of using the
+ // vec_accu_t alias for example) could possibly be frowned upon.
+
+ double x, y, z, length;
+
+ x = (double) in[0];
+ y = (double) in[1];
+ z = (double) in[2];
+
+ length = sqrt((x * x) + (y * y) + (z * z));
if (length == 0)
{
VectorClear (out);
return 0;
}
- ilength = 1.0f/length;
- out[0] = in[0]*ilength;
- out[1] = in[1]*ilength;
- out[2] = in[2]*ilength;
+ out[0] = (vec_t) (x / length);
+ out[1] = (vec_t) (y / length);
+ out[2] = (vec_t) (z / length);
- return length;
+ return (vec_t) length;
}
vec_t ColorNormalize( const vec3_t in, vec3_t out ) {
dst[i] = rot[i][0] * point[0] + rot[i][1] * point[1] + rot[i][2] * point[2];
}
}
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Below is double-precision math stuff. This was initially needed by the new
+// "base winding" code in q3map2 brush processing in order to fix the famous
+// "disappearing triangles" issue. These definitions can be used wherever extra
+// precision is needed.
+////////////////////////////////////////////////////////////////////////////////
+
+/*
+=================
+VectorLengthAccu
+=================
+*/
+vec_accu_t VectorLengthAccu(const vec3_accu_t v)
+{
+ return (vec_accu_t) sqrt((v[0] * v[0]) + (v[1] * v[1]) + (v[2] * v[2]));
+}
+
+/*
+=================
+DotProductAccu
+=================
+*/
+vec_accu_t DotProductAccu(const vec3_accu_t a, const vec3_accu_t b)
+{
+ return (a[0] * b[0]) + (a[1] * b[1]) + (a[2] * b[2]);
+}
+
+/*
+=================
+VectorSubtractAccu
+=================
+*/
+void VectorSubtractAccu(const vec3_accu_t a, const vec3_accu_t b, vec3_accu_t out)
+{
+ out[0] = a[0] - b[0];
+ out[1] = a[1] - b[1];
+ out[2] = a[2] - b[2];
+}
+
+/*
+=================
+VectorAddAccu
+=================
+*/
+void VectorAddAccu(const vec3_accu_t a, const vec3_accu_t b, vec3_accu_t out)
+{
+ out[0] = a[0] + b[0];
+ out[1] = a[1] + b[1];
+ out[2] = a[2] + b[2];
+}
+
+/*
+=================
+VectorCopyAccu
+=================
+*/
+void VectorCopyAccu(const vec3_accu_t in, vec3_accu_t out)
+{
+ out[0] = in[0];
+ out[1] = in[1];
+ out[2] = in[2];
+}
+
+/*
+=================
+VectorScaleAccu
+=================
+*/
+void VectorScaleAccu(const vec3_accu_t in, vec_accu_t scaleFactor, vec3_accu_t out)
+{
+ out[0] = in[0] * scaleFactor;
+ out[1] = in[1] * scaleFactor;
+ out[2] = in[2] * scaleFactor;
+}
+
+/*
+=================
+CrossProductAccu
+=================
+*/
+void CrossProductAccu(const vec3_accu_t a, const vec3_accu_t b, vec3_accu_t out)
+{
+ out[0] = (a[1] * b[2]) - (a[2] * b[1]);
+ out[1] = (a[2] * b[0]) - (a[0] * b[2]);
+ out[2] = (a[0] * b[1]) - (a[1] * b[0]);
+}
+
+/*
+=================
+Q_rintAccu
+=================
+*/
+vec_accu_t Q_rintAccu(vec_accu_t val)
+{
+ return (vec_accu_t) floor(val + 0.5);
+}
+
+/*
+=================
+VectorCopyAccuToRegular
+=================
+*/
+void VectorCopyAccuToRegular(const vec3_accu_t in, vec3_t out)
+{
+ out[0] = (vec_t) in[0];
+ out[1] = (vec_t) in[1];
+ out[2] = (vec_t) in[2];
+}
+
+/*
+=================
+VectorCopyRegularToAccu
+=================
+*/
+void VectorCopyRegularToAccu(const vec3_t in, vec3_accu_t out)
+{
+ out[0] = (vec_accu_t) in[0];
+ out[1] = (vec_accu_t) in[1];
+ out[2] = (vec_accu_t) in[2];
+}
+
+/*
+=================
+VectorNormalizeAccu
+=================
+*/
+vec_accu_t VectorNormalizeAccu(const vec3_accu_t in, vec3_accu_t out)
+{
+ // The sqrt() function takes double as an input and returns double as an
+ // output according the the man pages on Debian and on FreeBSD. Therefore,
+ // I don't see a reason why using a double outright (instead of using the
+ // vec_accu_t alias for example) could possibly be frowned upon.
+
+ vec_accu_t length;
+
+ length = (vec_accu_t) sqrt((in[0] * in[0]) + (in[1] * in[1]) + (in[2] * in[2]));
+ if (length == 0)
+ {
+ VectorClear(out);
+ return 0;
+ }
+
+ out[0] = in[0] / length;
+ out[1] = in[1] / length;
+ out[2] = in[2] / length;
+
+ return length;
+}
return w;
}
+/*
+=============
+AllocWindingAccu
+=============
+*/
+winding_accu_t *AllocWindingAccu(int points)
+{
+ winding_accu_t *w;
+ int s;
+
+ if (points >= MAX_POINTS_ON_WINDING)
+ Error("AllocWindingAccu failed: MAX_POINTS_ON_WINDING exceeded");
+
+ if (numthreads == 1)
+ {
+ // At the time of this writing, these statistics were not used in any way.
+ c_winding_allocs++;
+ c_winding_points += points;
+ c_active_windings++;
+ if (c_active_windings > c_peak_windings)
+ c_peak_windings = c_active_windings;
+ }
+ s = sizeof(vec_accu_t) * 3 * points + sizeof(int);
+ w = safe_malloc(s);
+ memset(w, 0, s);
+ return w;
+}
+
+/*
+=============
+FreeWinding
+=============
+*/
void FreeWinding (winding_t *w)
{
+ if (!w) Error("FreeWinding: winding is NULL");
+
if (*(unsigned *)w == 0xdeaddead)
Error ("FreeWinding: freed a freed winding");
*(unsigned *)w = 0xdeaddead;
free (w);
}
+/*
+=============
+FreeWindingAccu
+=============
+*/
+void FreeWindingAccu(winding_accu_t *w)
+{
+ if (!w) Error("FreeWindingAccu: winding is NULL");
+
+ if (*((unsigned *) w) == 0xdeaddead)
+ Error("FreeWindingAccu: freed a freed winding");
+ *((unsigned *) w) = 0xdeaddead;
+
+ if (numthreads == 1)
+ c_active_windings--;
+ free(w);
+}
+
/*
============
RemoveColinearPoints
VectorScale (center, scale, center);
}
+/*
+=================
+BaseWindingForPlaneAccu
+=================
+*/
+winding_accu_t *BaseWindingForPlaneAccu(vec3_t normal, vec_t dist)
+{
+ // The goal in this function is to replicate the behavior of the original BaseWindingForPlane()
+ // function (see below) but at the same time increasing accuracy substantially.
+
+ // The original code gave a preference for the vup vector to start out as (0, 0, 1), unless the
+ // normal had a dominant Z value, in which case vup started out as (1, 0, 0). After that, vup
+ // was "bent" [along the plane defined by normal and vup] to become perpendicular to normal.
+ // After that the vright vector was computed as the cross product of vup and normal.
+
+ // I'm constructing the winding polygon points in a fashion similar to the method used in the
+ // original function. Orientation is the same. The size of the winding polygon, however, is
+ // variable in this function (depending on the angle of normal), and is larger (by about a factor
+ // of 2) than the winding polygon in the original function.
+
+ int x, i;
+ vec_t max, v;
+ vec3_accu_t vright, vup, org, normalAccu;
+ winding_accu_t *w;
+
+ // One of the components of normal must have a magnitiude greater than this value,
+ // otherwise normal is not a unit vector. This is a little bit of inexpensive
+ // partial error checking we can do.
+ max = 0.56; // 1 / sqrt(1^2 + 1^2 + 1^2) = 0.577350269
+
+ x = -1;
+ for (i = 0; i < 3; i++) {
+ v = (vec_t) fabs(normal[i]);
+ if (v > max) {
+ x = i;
+ max = v;
+ }
+ }
+ if (x == -1) Error("BaseWindingForPlaneAccu: no dominant axis found because normal is too short");
+
+ switch (x) {
+ case 0: // Fall through to next case.
+ case 1:
+ vright[0] = (vec_accu_t) -normal[1];
+ vright[1] = (vec_accu_t) normal[0];
+ vright[2] = 0;
+ break;
+ case 2:
+ vright[0] = 0;
+ vright[1] = (vec_accu_t) -normal[2];
+ vright[2] = (vec_accu_t) normal[1];
+ break;
+ }
+
+ // vright and normal are now perpendicular; you can prove this by taking their
+ // dot product and seeing that it's always exactly 0 (with no error).
+
+ // NOTE: vright is NOT a unit vector at this point. vright will have length
+ // not exceeding 1.0. The minimum length that vright can achieve happens when,
+ // for example, the Z and X components of the normal input vector are equal,
+ // and when normal's Y component is zero. In that case Z and X of the normal
+ // vector are both approximately 0.70711. The resulting vright vector in this
+ // case will have a length of 0.70711.
+
+ // We're relying on the fact that MAX_WORLD_COORD is a power of 2 to keep
+ // our calculation precise and relatively free of floating point error.
+ // [However, the code will still work fine if that's not the case.]
+ VectorScaleAccu(vright, ((vec_accu_t) MAX_WORLD_COORD) * 4.0, vright);
+
+ // At time time of this writing, MAX_WORLD_COORD was 65536 (2^16). Therefore
+ // the length of vright at this point is at least 185364. In comparison, a
+ // corner of the world at location (65536, 65536, 65536) is distance 113512
+ // away from the origin.
+
+ VectorCopyRegularToAccu(normal, normalAccu);
+ CrossProductAccu(normalAccu, vright, vup);
+
+ // vup now has length equal to that of vright.
+
+ VectorScaleAccu(normalAccu, (vec_accu_t) dist, org);
+
+ // org is now a point on the plane defined by normal and dist. Furthermore,
+ // org, vright, and vup are pairwise perpendicular. Now, the 4 vectors
+ // { (+-)vright + (+-)vup } have length that is at least sqrt(185364^2 + 185364^2),
+ // which is about 262144. That length lies outside the world, since the furthest
+ // point in the world has distance 113512 from the origin as mentioned above.
+ // Also, these 4 vectors are perpendicular to the org vector. So adding them
+ // to org will only increase their length. Therefore the 4 points defined below
+ // all lie outside of the world. Furthermore, it can be easily seen that the
+ // edges connecting these 4 points (in the winding_accu_t below) lie completely
+ // outside the world. sqrt(262144^2 + 262144^2)/2 = 185363, which is greater than
+ // 113512.
+
+ w = AllocWindingAccu(4);
+
+ VectorSubtractAccu(org, vright, w->p[0]);
+ VectorAddAccu(w->p[0], vup, w->p[0]);
+
+ VectorAddAccu(org, vright, w->p[1]);
+ VectorAddAccu(w->p[1], vup, w->p[1]);
+
+ VectorAddAccu(org, vright, w->p[2]);
+ VectorSubtractAccu(w->p[2], vup, w->p[2]);
+
+ VectorSubtractAccu(org, vright, w->p[3]);
+ VectorSubtractAccu(w->p[3], vup, w->p[3]);
+
+ w->numpoints = 4;
+
+ return w;
+}
+
/*
=================
BaseWindingForPlane
+
+Original BaseWindingForPlane() function that has serious accuracy problems. Here is why.
+The base winding is computed as a rectangle with very large coordinates. These coordinates
+are in the range 2^17 or 2^18. "Epsilon" (meaning the distance between two adjacent numbers)
+at these magnitudes in 32 bit floating point world is about 0.02. So for example, if things
+go badly (by bad luck), then the whole plane could be shifted by 0.02 units (its distance could
+be off by that much). Then if we were to compute the winding of this plane and another of
+the brush's planes met this winding at a very acute angle, that error could multiply to around
+0.1 or more when computing the final vertex coordinates of the winding. 0.1 is a very large
+error, and can lead to all sorts of disappearing triangle problems.
=================
*/
winding_t *BaseWindingForPlane (vec3_t normal, vec_t dist)
int size;
winding_t *c;
+ if (!w) Error("CopyWinding: winding is NULL");
+
c = AllocWinding (w->numpoints);
size = (int)((size_t)((winding_t *)0)->p[w->numpoints]);
memcpy (c, w, size);
return c;
}
+/*
+==================
+CopyWindingAccuIncreaseSizeAndFreeOld
+==================
+*/
+winding_accu_t *CopyWindingAccuIncreaseSizeAndFreeOld(winding_accu_t *w)
+{
+ int i;
+ winding_accu_t *c;
+
+ if (!w) Error("CopyWindingAccuIncreaseSizeAndFreeOld: winding is NULL");
+
+ c = AllocWindingAccu(w->numpoints + 1);
+ c->numpoints = w->numpoints;
+ for (i = 0; i < c->numpoints; i++)
+ {
+ VectorCopyAccu(w->p[i], c->p[i]);
+ }
+ FreeWindingAccu(w);
+ return c;
+}
+
+/*
+==================
+CopyWindingAccuToRegular
+==================
+*/
+winding_t *CopyWindingAccuToRegular(winding_accu_t *w)
+{
+ int i;
+ winding_t *c;
+
+ if (!w) Error("CopyWindingAccuToRegular: winding is NULL");
+
+ c = AllocWinding(w->numpoints);
+ c->numpoints = w->numpoints;
+ for (i = 0; i < c->numpoints; i++)
+ {
+ VectorCopyAccuToRegular(w->p[i], c->p[i]);
+ }
+ return c;
+}
+
/*
==================
ReverseWinding
}
+/*
+=============
+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, normalAccu;
+
+ // 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 are 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.
+
+ // Here is a further discussion regarding the choice of smallestEpsilonAllowed.
+ // In the 32 float world (we can assume vec_t is that), the "epsilon around 1.0" is
+ // 0.00000011921. In the 64 bit float world (we can assume vec_accu_t is that), the
+ // "epsilon around 1.0" is 0.00000000000000022204. (By the way these two epsilons
+ // are defined as VEC_SMALLEST_EPSILON_AROUND_ONE VEC_ACCU_SMALLEST_EPSILON_AROUND_ONE
+ // respectively.) If you divide the first by the second, you get approximately
+ // 536,885,246. Dividing that number by 200,000 (a typical base winding coordinate)
+ // gives 2684. So in other words, if our smallestEpsilonAllowed was chosen as exactly
+ // VEC_SMALLEST_EPSILON_AROUND_ONE, you would be guaranteed at least 2000 "ticks" in
+ // 64-bit land inside of the epsilon for all numbers we're dealing with.
+
+ static const vec_accu_t smallestEpsilonAllowed = ((vec_accu_t) VEC_SMALLEST_EPSILON_AROUND_ONE) * 0.5;
+ if (crudeEpsilon < smallestEpsilonAllowed) fineEpsilon = smallestEpsilonAllowed;
+ else fineEpsilon = (vec_accu_t) crudeEpsilon;
+
+ in = *inout;
+ counts[0] = counts[1] = counts[2] = 0;
+ VectorCopyRegularToAccu(normal, normalAccu);
+
+ for (i = 0; i < in->numpoints; i++)
+ {
+ dists[i] = DotProductAccu(in->p[i], normalAccu) - 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];
+
+ // I'm wondering if whatever code that handles duplicate planes is robust enough
+ // that we never get a case where two nearly equal planes result in 2 NULL windings
+ // due to the 'if' statement below. TODO: Investigate this.
+ if (!counts[SIDE_FRONT]) {
+ FreeWindingAccu(in);
+ *inout = NULL;
+ return;
+ }
+ if (!counts[SIDE_BACK]) {
+ return; // Winding is unmodified.
+ }
+
+ // NOTE: The least number of points that a winding can have at this point is 2.
+ // In that case, one point is SIDE_FRONT and the other is SIDE_BACK.
+
+ maxpts = counts[SIDE_FRONT] + 2; // We dynamically expand if this is too small.
+ f = AllocWindingAccu(maxpts);
+
+ for (i = 0; i < in->numpoints; i++)
+ {
+ p1 = in->p[i];
+
+ if (sides[i] == SIDE_ON || sides[i] == SIDE_FRONT)
+ {
+ if (f->numpoints >= MAX_POINTS_ON_WINDING)
+ Error("ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING");
+ if (f->numpoints >= maxpts) // This will probably never happen.
+ {
+ Sys_FPrintf(SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n");
+ f = CopyWindingAccuIncreaseSizeAndFreeOld(f);
+ maxpts++;
+ }
+ VectorCopyAccu(p1, f->p[f->numpoints]);
+ f->numpoints++;
+ if (sides[i] == SIDE_ON) continue;
+ }
+ 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 >= MAX_POINTS_ON_WINDING)
+ Error("ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING");
+ if (f->numpoints >= maxpts) // This will probably never happen.
+ {
+ Sys_FPrintf(SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n");
+ f = CopyWindingAccuIncreaseSizeAndFreeOld(f);
+ maxpts++;
+ }
+ VectorCopyAccu(mid, f->p[f->numpoints]);
+ f->numpoints++;
+ }
+
+ FreeWindingAccu(in);
+ *inout = f;
+}
+
/*
=============
ChopWindingInPlace
// frees the original if clipped
void pw(winding_t *w);
+
+
+///////////////////////////////////////////////////////////////////////////////////////
+// Below is double-precision stuff. This was initially needed by the base winding code
+// in q3map2 brush processing.
+///////////////////////////////////////////////////////////////////////////////////////
+
+typedef struct
+{
+ int numpoints;
+ vec3_accu_t p[4]; // variable sized
+} 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);
+winding_t *CopyWindingAccuToRegular(winding_accu_t *w);
+void FreeWindingAccu(winding_accu_t *w);
}
}
+/*
+==================
+SnapWeldVectorAccu
+
+Welds two vectors into a third, taking into account nearest-to-integer
+instead of averaging.
+==================
+*/
+void SnapWeldVectorAccu(vec3_accu_t a, vec3_accu_t b, vec3_accu_t out)
+{
+ // I'm just preserving what I think was the intended logic of the original
+ // SnapWeldVector(). I'm not actually sure where this function should even
+ // be used. I'd like to know which kinds of problems this function addresses.
+
+ // TODO: I thought we're snapping all coordinates to nearest 1/8 unit?
+ // So what is natural about snapping to the nearest integer? Maybe we should
+ // be snapping to the nearest 1/8 unit instead?
+
+ int i;
+ vec_accu_t ai, bi, ad, bd;
+
+ if (a == NULL || b == NULL || out == NULL)
+ Error("SnapWeldVectorAccu: NULL argument");
+
+ for (i = 0; i < 3; i++)
+ {
+ ai = Q_rintAccu(a[i]);
+ bi = Q_rintAccu(b[i]);
+ ad = fabs(ai - a[i]);
+ bd = fabs(bi - b[i]);
+
+ if (ad < bd)
+ {
+ if (ad < SNAP_EPSILON) out[i] = ai;
+ else out[i] = a[i];
+ }
+ else
+ {
+ if (bd < SNAP_EPSILON) out[i] = bi;
+ else out[i] = b[i];
+ }
+ }
+}
+
/*
return valid;
}
+/*
+==================
+FixWindingAccu
+
+Removes degenerate edges (edges that are too short) from a winding.
+Returns qtrue if the winding has been altered by this function.
+Returns qfalse if the winding is untouched by this function.
+
+It's advised that you check the winding after this function exits to make
+sure it still has at least 3 points. If that is not the case, the winding
+cannot be considered valid. The winding may degenerate to one or two points
+if the some of the winding's points are close together.
+==================
+*/
+qboolean FixWindingAccu(winding_accu_t *w)
+{
+ int i, j, k;
+ vec3_accu_t vec;
+ vec_accu_t dist;
+ qboolean done, altered;
+
+ if (w == NULL) Error("FixWindingAccu: NULL argument");
+ altered = qfalse;
+ while (qtrue)
+ {
+ if (w->numpoints < 2) break; // Don't remove the only remaining point.
+ done = qtrue;
+ for (i = 0; i < w->numpoints; i++)
+ {
+ j = (((i + 1) == w->numpoints) ? 0 : (i + 1));
+ VectorSubtractAccu(w->p[i], w->p[j], vec);
+ dist = VectorLengthAccu(vec);
+ if (dist < DEGENERATE_EPSILON)
+ {
+ // TODO: I think the "snap weld vector" was written before
+ // some of the math precision fixes, and its purpose was
+ // probably to address math accuracy issues. We can think
+ // about changing the logic here. Maybe once plane distance
+ // gets 64 bits, we can look at it then.
+ SnapWeldVectorAccu(w->p[i], w->p[j], vec);
+ VectorCopyAccu(vec, w->p[i]);
+ for (k = j + 1; k < w->numpoints; k++)
+ {
+ VectorCopyAccu(w->p[k], w->p[k - 1]);
+ }
+ w->numpoints--;
+ altered = qtrue;
+ // The only way to finish off fixing the winding consistently and
+ // accurately is by fixing the winding all over again. For example,
+ // the point at index i and the point at index i-1 could now be
+ // less than the epsilon distance apart. There are too many special
+ // case problems we'd need to handle if we didn't start from the
+ // beginning.
+ done = qfalse;
+ break; // This will cause us to return to the "while" loop.
+ }
+ }
+ if (done) break;
+ }
+ return altered;
+}
/*
qboolean CreateBrushWindings( brush_t *brush )
{
int i, j;
+#if EXPERIMENTAL_HIGH_PRECISION_MATH_Q3MAP2_FIXES
+ winding_accu_t *w;
+#else
winding_t *w;
+#endif
side_t *side;
plane_t *plane;
plane = &mapplanes[ side->planenum ];
/* make huge winding */
+#if EXPERIMENTAL_HIGH_PRECISION_MATH_Q3MAP2_FIXES
+ w = BaseWindingForPlaneAccu(plane->normal, plane->dist);
+#else
w = BaseWindingForPlane( plane->normal, plane->dist );
+#endif
/* walk the list of brush sides */
for( j = 0; j < brush->numsides && w != NULL; j++ )
if( brush->sides[ j ].bevel )
continue;
plane = &mapplanes[ brush->sides[ j ].planenum ^ 1 ];
+#if EXPERIMENTAL_HIGH_PRECISION_MATH_Q3MAP2_FIXES
+ ChopWindingInPlaceAccu(&w, plane->normal, plane->dist, 0);
+#else
ChopWindingInPlace( &w, plane->normal, plane->dist, 0 ); // CLIP_EPSILON );
+#endif
/* ydnar: fix broken windings that would generate trifans */
+#if EXPERIMENTAL_HIGH_PRECISION_MATH_Q3MAP2_FIXES
+ // I think it's better to FixWindingAccu() once after we chop with all planes
+ // so that error isn't multiplied. There is nothing natural about welding
+ // the points unless they are the final endpoints. ChopWindingInPlaceAccu()
+ // is able to handle all kinds of degenerate windings.
+#else
FixWinding( w );
+#endif
}
/* set side winding */
+#if EXPERIMENTAL_HIGH_PRECISION_MATH_Q3MAP2_FIXES
+ if (w != NULL)
+ {
+ FixWindingAccu(w);
+ if (w->numpoints < 3)
+ {
+ FreeWindingAccu(w);
+ w = NULL;
+ }
+ }
+ side->winding = (w ? CopyWindingAccuToRegular(w) : NULL);
+ if (w) FreeWindingAccu(w);
+#else
side->winding = w;
+#endif
}
/* find brush bounds */
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]);
ydnar: replaced with variable epsilon for djbob
*/
-#define NORMAL_EPSILON 0.00001
-#define DIST_EPSILON 0.01
-
qboolean PlaneEqual( plane_t *p, vec3_t normal, vec_t dist )
{
float ne, de;
de = distanceEpsilon;
/* compare */
- if( fabs( p->dist - dist ) <= de &&
- fabs( p->normal[ 0 ] - normal[ 0 ] ) <= ne &&
- fabs( p->normal[ 1 ] - normal[ 1 ] ) <= ne &&
- fabs( p->normal[ 2 ] - normal[ 2 ] ) <= ne )
+ // We check equality of each component since we're using '<', not '<='
+ // (the epsilons may be zero). We want to use '<' intead of '<=' to be
+ // consistent with the true meaning of "epsilon", and also because other
+ // parts of the code uses this inequality.
+ if ((p->dist == dist || fabs(p->dist - dist) < de) &&
+ (p->normal[0] == normal[0] || fabs(p->normal[0] - normal[0]) < ne) &&
+ (p->normal[1] == normal[1] || fabs(p->normal[1] - normal[1]) < ne) &&
+ (p->normal[2] == normal[2] || fabs(p->normal[2] - normal[2]) < ne))
return qtrue;
/* different */
/*
SnapNormal()
-snaps a near-axial normal vector
+Snaps a near-axial normal vector.
+Returns qtrue if and only if the normal was adjusted.
*/
-void SnapNormal( vec3_t normal )
+qboolean SnapNormal( vec3_t normal )
{
+#if EXPERIMENTAL_SNAP_NORMAL_FIX
+ int i;
+ qboolean adjusted = qfalse;
+
+ // A change from the original SnapNormal() is that we snap each
+ // component that's close to 0. So for example if a normal is
+ // (0.707, 0.707, 0.0000001), it will get snapped to lie perfectly in the
+ // XY plane (its Z component will be set to 0 and its length will be
+ // normalized). The original SnapNormal() didn't snap such vectors - it
+ // only snapped vectors that were near a perfect axis.
+
+ for (i = 0; i < 3; i++)
+ {
+ if (normal[i] != 0.0 && -normalEpsilon < normal[i] && normal[i] < normalEpsilon)
+ {
+ normal[i] = 0.0;
+ adjusted = qtrue;
+ }
+ }
+
+ if (adjusted)
+ {
+ VectorNormalize(normal, normal);
+ return qtrue;
+ }
+ return qfalse;
+#else
int i;
+ // I would suggest that you uncomment the following code and look at the
+ // results:
+
+ /*
+ Sys_Printf("normalEpsilon is %f\n", normalEpsilon);
+ for (i = 0;; i++)
+ {
+ normal[0] = 1.0;
+ normal[1] = 0.0;
+ normal[2] = i * 0.000001;
+ VectorNormalize(normal, normal);
+ if (1.0 - normal[0] >= normalEpsilon) {
+ Sys_Printf("(%f %f %f)\n", normal[0], normal[1], normal[2]);
+ Error("SnapNormal: test completed");
+ }
+ }
+ */
+
+ // When the normalEpsilon is 0.00001, the loop will break out when normal is
+ // (0.999990 0.000000 0.004469). In other words, this is the vector closest
+ // to axial that will NOT be snapped. Anything closer will be snaped. Now,
+ // 0.004469 is close to 1/225. The length of a circular quarter-arc of radius
+ // 1 is PI/2, or about 1.57. And 0.004469/1.57 is about 0.0028, or about
+ // 1/350. Expressed a different way, 1/350 is also about 0.26/90.
+ // This means is that a normal with an angle that is within 1/4 of a degree
+ // from axial will be "snapped". My belief is that the person who wrote the
+ // code below did not intend it this way. I think the person intended that
+ // the epsilon be measured against the vector components close to 0, not 1.0.
+ // I think the logic should be: if 2 of the normal components are within
+ // epsilon of 0, then the vector can be snapped to be perfectly axial.
+ // We may consider adjusting the epsilon to a larger value when we make this
+ // code fix.
+
for( i = 0; i < 3; i++ )
{
if( fabs( normal[ i ] - 1 ) < normalEpsilon )
{
VectorClear( normal );
normal[ i ] = 1;
- break;
+ return qtrue;
}
if( fabs( normal[ i ] - -1 ) < normalEpsilon )
{
VectorClear( normal );
normal[ i ] = -1;
- break;
+ return qtrue;
}
}
+ return qfalse;
+#endif
}
*/
SnapNormal( normal );
+ // TODO: Rambetter has some serious comments here as well. First off,
+ // in the case where a normal is non-axial, there is nothing special
+ // about integer distances. I would think that snapping a distance might
+ // make sense for axial normals, but I'm not so sure about snapping
+ // non-axial normals. A shift by 0.01 in a plane, multiplied by a clipping
+ // against another plane that is 5 degrees off, and we introduce 0.1 error
+ // easily. A 0.1 error in a vertex is where problems start to happen, such
+ // as disappearing triangles.
+
+ // Second, assuming we have snapped the normal above, let's say that the
+ // plane we just snapped was defined for some points that are actually
+ // quite far away from normal * dist. Well, snapping the normal in this
+ // case means that we've just moved those points by potentially many units!
+ // Therefore, if we are going to snap the normal, we need to know the
+ // points we're snapping for so that the plane snaps with those points in
+ // mind (points remain close to the plane).
+
+ // I would like to know exactly which problems SnapPlane() is trying to
+ // solve so that we can better engineer it (I'm not saying that SnapPlane()
+ // should be removed altogether). Fix all this snapping code at some point!
+
if( fabs( *dist - Q_rint( *dist ) ) < distanceEpsilon )
*dist = Q_rint( *dist );
}
+/*
+SnapPlaneImproved()
+snaps a plane to normal/distance epsilons, improved code
+*/
+void SnapPlaneImproved(vec3_t normal, vec_t *dist, int numPoints, const vec3_t *points)
+{
+ int i;
+ vec3_t center;
+ vec_t distNearestInt;
+
+ if (SnapNormal(normal))
+ {
+ if (numPoints > 0)
+ {
+ // Adjust the dist so that the provided points don't drift away.
+ VectorClear(center);
+ for (i = 0; i < numPoints; i++)
+ {
+ VectorAdd(center, points[i], center);
+ }
+ for (i = 0; i < 3; i++) { center[i] = center[i] / numPoints; }
+ *dist = DotProduct(normal, center);
+ }
+ }
+
+ if (VectorIsOnAxis(normal))
+ {
+ // Only snap distance if the normal is an axis. Otherwise there
+ // is nothing "natural" about snapping the distance to an integer.
+ distNearestInt = Q_rint(*dist);
+ if (-distanceEpsilon < *dist - distNearestInt && *dist - distNearestInt < distanceEpsilon)
+ {
+ *dist = distNearestInt;
+ }
+ }
+}
/*
vec_t d;
- /* hash the plane */
+#if EXPERIMENTAL_SNAP_PLANE_FIX
+ SnapPlaneImproved(normal, &dist, numPoints, (const vec3_t *) points);
+#else
SnapPlane( normal, &dist );
+#endif
+ /* hash the plane */
hash = (PLANE_HASHES - 1) & (int) fabs( dist );
/* search the border bins as well */
/* ydnar: test supplied points against this plane */
for( j = 0; j < numPoints; j++ )
{
+ // NOTE: When dist approaches 2^16, the resolution of 32 bit floating
+ // point number is greatly decreased. The distanceEpsilon cannot be
+ // very small when world coordinates extend to 2^16. Making the
+ // dot product here in 64 bit land will not really help the situation
+ // because the error will already be carried in dist.
d = DotProduct( points[ j ], normal ) - dist;
- if( fabs( d ) > distanceEpsilon )
- break;
+ d = fabs(d);
+ if (d != 0.0 && d >= distanceEpsilon)
+ break; // Point is too far from plane.
}
/* found a matching plane */
int i;
plane_t *p;
-
+#if EXPERIMENTAL_SNAP_PLANE_FIX
+ SnapPlaneImproved(normal, &dist, numPoints, (const vec3_t *) points);
+#else
SnapPlane( normal, &dist );
+#endif
for( i = 0, p = mapplanes; i < nummapplanes; i++, p++ )
{
if( PlaneEqual( p, normal, dist ) )
return i;
+ // TODO: Note that the non-USE_HASHING code does not compute epsilons
+ // for the provided points. It should do that. I think this code
+ // is unmaintained because nobody sets USE_HASHING to off.
}
return CreateNewFloatPlane( normal, dist );
int MapPlaneFromPoints( vec3_t *p )
{
+#if EXPERIMENTAL_HIGH_PRECISION_MATH_Q3MAP2_FIXES
+ vec3_accu_t paccu[3];
+ vec3_accu_t t1, t2, normalAccu;
+ vec3_t normal;
+ vec_t dist;
+
+ VectorCopyRegularToAccu(p[0], paccu[0]);
+ VectorCopyRegularToAccu(p[1], paccu[1]);
+ VectorCopyRegularToAccu(p[2], paccu[2]);
+
+ VectorSubtractAccu(paccu[0], paccu[1], t1);
+ VectorSubtractAccu(paccu[2], paccu[1], t2);
+ CrossProductAccu(t1, t2, normalAccu);
+ VectorNormalizeAccu(normalAccu, normalAccu);
+ // TODO: A 32 bit float for the plane distance isn't enough resolution
+ // if the plane is 2^16 units away from the origin (the "epsilon" approaches
+ // 0.01 in that case).
+ dist = (vec_t) DotProductAccu(paccu[0], normalAccu);
+ VectorCopyAccuToRegular(normalAccu, normal);
+
+ return FindFloatPlane(normal, dist, 3, p);
+#else
vec3_t t1, t2, normal;
vec_t dist;
/* store the plane */
return FindFloatPlane( normal, dist, 3, p );
+#endif
}
------------------------------------------------------------------------------- */
+/* temporary hacks and tests (please keep off in SVN to prevent anyone's legacy map from screwing up) */
+#define EXPERIMENTAL_HIGH_PRECISION_MATH_Q3MAP2_FIXES 0
+#define EXPERIMENTAL_SNAP_NORMAL_FIX 0
+#define EXPERIMENTAL_SNAP_PLANE_FIX 0
+
/* general */
#define MAX_QPATH 64
Q_EXTERN qboolean debugInset Q_ASSIGN( qfalse );
Q_EXTERN qboolean debugPortals Q_ASSIGN( qfalse );
+#if EXPERIMENTAL_SNAP_NORMAL_FIX
+// Increasing the normalEpsilon to compensate for new logic in SnapNormal(), where
+// this epsilon is now used to compare against 0 components instead of the 1 or -1
+// components. Unfortunately, normalEpsilon is also used in PlaneEqual(). So changing
+// this will affect anything that calls PlaneEqual() as well (which are, at the time
+// of this writing, FindFloatPlane() and AddBrushBevels()).
+Q_EXTERN double normalEpsilon Q_ASSIGN(0.00005);
+#else
Q_EXTERN double normalEpsilon Q_ASSIGN( 0.00001 );
+#endif
+
+#if EXPERIMENTAL_HIGH_PRECISION_MATH_Q3MAP2_FIXES
+// NOTE: This distanceEpsilon is too small if parts of the map are at maximum world
+// extents (in the range of plus or minus 2^16). The smallest epsilon at values
+// close to 2^16 is about 0.007, which is greater than distanceEpsilon. Therefore,
+// maps should be constrained to about 2^15, otherwise slightly undesirable effects
+// may result. The 0.01 distanceEpsilon used previously is just too coarse in my
+// opinion. The real fix for this problem is to have 64 bit distances and then make
+// this epsilon even smaller, or to constrain world coordinates to plus minus 2^15
+// (or even 2^14).
+Q_EXTERN double distanceEpsilon Q_ASSIGN(0.005);
+#else
Q_EXTERN double distanceEpsilon Q_ASSIGN( 0.01 );
+#endif
/* bsp */
projects / xonotic / netradiant.git / commitdiff