From: rambetter Date: Mon, 10 Jan 2011 06:15:47 +0000 (+0000) Subject: Importing code changes for q3map2 from Rambetter-math-fix-experiments branch X-Git-Url: https://git.rm.cloudns.org/?a=commitdiff_plain;h=7006e9eead128058dfe0a5efd121fe3e41e71fc1;p=xonotic%2Fnetradiant.git Importing code changes for q3map2 from Rambetter-math-fix-experiments branch into trunk. Right now all the new code that fixes problems is turned off. There are three new #defines in q3map2.h: EXPERIMENTAL_HIGH_PRECISION_MATH_Q3MAP2_FIXES, EXPERIMENTAL_SNAP_NORMAL_FIX, and EXPERIMENTAL_SNAP_PLANE_FIX. All of these are currently set to 0, which means don't enable that new code. You can easily edit these to be 1 in order to enable the new code. There are very very minor changes to the code even with these three #defines disabled. They are as follows. - In PlaneEqual() in map.c, now considering deltas equal to given epsilon values as "far enough to be different". Previously, the '<=' operation was used, now '<' is being used. - In FindFloatPlane() in map.c, considering delta equal to distanceEpsilon (for plane distance) to be sufficiently far away. Before, delta had to be strictly greater than distanceEpsilon. - VectorNormalize() in mathlib.c is more accurate now. This change itself causes at least one regression test to succeed. The previous implementation of VectorNormalize() caused excessive errors to be introduced due to sloppy arithmetic. Note, the epsilon changes account for the possibility that the epsilons are set to 0.0 on the command-line. git-svn-id: https://zerowing.idsoftware.com/svn/radiant/GtkRadiant/trunk@416 8a3a26a2-13c4-0310-b231-cf6edde360e5 --- diff --git a/libs/mathlib.h b/libs/mathlib.h index 44bf01c9..cac3597b 100644 --- a/libs/mathlib.h +++ b/libs/mathlib.h @@ -24,6 +24,7 @@ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA // mathlib.h #include +#include #include "bytebool.h" @@ -37,6 +38,12 @@ typedef vec_t vec3_t[3]; 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 @@ -75,6 +82,9 @@ qboolean VectorCompare (vec3_t v1, vec3_t v2); #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 ); @@ -298,6 +308,50 @@ vec_t ray_intersect_point(const ray_t *ray, const vec3_t point, vec_t epsilon, v /*! 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 diff --git a/libs/mathlib/mathlib.c b/libs/mathlib/mathlib.c index b9b7b869..0738c9c4 100644 --- a/libs/mathlib/mathlib.c +++ b/libs/mathlib/mathlib.c @@ -26,6 +26,54 @@ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 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 @@ -127,21 +175,30 @@ void _VectorCopy (vec3_t in, vec3_t out) } 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 ) { @@ -591,3 +648,153 @@ void RotatePointAroundVector( vec3_t dst, const vec3_t dir, const vec3_t point, 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; +} diff --git a/tools/quake3/common/polylib.c b/tools/quake3/common/polylib.c index 75caec6b..78444f94 100644 --- a/tools/quake3/common/polylib.c +++ b/tools/quake3/common/polylib.c @@ -73,8 +73,43 @@ winding_t *AllocWinding (int points) 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; @@ -84,6 +119,24 @@ void FreeWinding (winding_t *w) 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 @@ -201,9 +254,131 @@ void WindingCenter (winding_t *w, vec3_t center) 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) @@ -285,12 +460,57 @@ winding_t *CopyWinding (winding_t *w) 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 @@ -424,6 +644,147 @@ 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, 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 diff --git a/tools/quake3/common/polylib.h b/tools/quake3/common/polylib.h index 58b21d5a..5e560d9d 100644 --- a/tools/quake3/common/polylib.h +++ b/tools/quake3/common/polylib.h @@ -55,3 +55,20 @@ void ChopWindingInPlace (winding_t **w, vec3_t normal, vec_t dist, vec_t epsilon // 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); diff --git a/tools/quake3/q3map2/brush.c b/tools/quake3/q3map2/brush.c index 2517ca63..12754085 100644 --- a/tools/quake3/q3map2/brush.c +++ b/tools/quake3/q3map2/brush.c @@ -277,6 +277,50 @@ void SnapWeldVector( vec3_t a, vec3_t b, vec3_t out ) } } +/* +================== +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]; + } + } +} + /* @@ -338,10 +382,71 @@ qboolean FixWinding( winding_t *w ) 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; +} /* @@ -353,7 +458,11 @@ returns false if the brush doesn't enclose a valid volume 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; @@ -366,7 +475,11 @@ qboolean CreateBrushWindings( brush_t *brush ) 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++ ) @@ -378,14 +491,39 @@ qboolean CreateBrushWindings( brush_t *brush ) 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 */ @@ -506,6 +644,8 @@ void WriteBSPBrushMap( char *name, brush_t *list ) fprintf (f, "{\n"); for (i=0,s=list->sides ; inumsides ; 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]); diff --git a/tools/quake3/q3map2/map.c b/tools/quake3/q3map2/map.c index 1e2f3267..330e2de5 100644 --- a/tools/quake3/q3map2/map.c +++ b/tools/quake3/q3map2/map.c @@ -59,9 +59,6 @@ PlaneEqual() 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; @@ -72,10 +69,14 @@ qboolean PlaneEqual( plane_t *p, vec3_t normal, vec_t dist ) 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 */ @@ -153,28 +154,91 @@ int CreateNewFloatPlane (vec3_t normal, vec_t dist) /* 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 } @@ -195,10 +259,67 @@ void SnapPlane( vec3_t normal, vec_t *dist ) */ 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; + } + } +} /* @@ -217,8 +338,12 @@ int FindFloatPlane( vec3_t normal, vec_t dist, int numPoints, vec3_t *points ) 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 */ @@ -237,9 +362,15 @@ int FindFloatPlane( vec3_t normal, vec_t dist, int numPoints, vec3_t *points ) /* 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 */ @@ -258,12 +389,18 @@ int FindFloatPlane( vec3_t normal, vec_t dist, int numPoints, vec3_t *points ) 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 ); @@ -280,6 +417,28 @@ takes 3 points and finds the plane they lie in 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; @@ -295,6 +454,7 @@ int MapPlaneFromPoints( vec3_t *p ) /* store the plane */ return FindFloatPlane( normal, dist, 3, p ); +#endif } diff --git a/tools/quake3/q3map2/q3map2.h b/tools/quake3/q3map2/q3map2.h index 9145037b..e3de7b7a 100644 --- a/tools/quake3/q3map2/q3map2.h +++ b/tools/quake3/q3map2/q3map2.h @@ -122,6 +122,11 @@ constants ------------------------------------------------------------------------------- */ +/* 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 @@ -1925,8 +1930,30 @@ Q_EXTERN qboolean debugSurfaces Q_ASSIGN( qfalse ); 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 */