From 94cb905da9c38e65cff6dcd9a56874c33c367df2 Mon Sep 17 00:00:00 2001 From: Rudolf Polzer Date: Tue, 11 Jan 2011 14:39:55 +0100 Subject: [PATCH] ::zerowing-base=422 --- libs/mathlib.h | 54 +++ libs/mathlib/mathlib.c | 221 ++++++++++- .../disappearing_sliver/winding_logging.patch | 54 ++- .../q3map2/disappearing_sliver3/NOTES.txt | 36 ++ tools/quake3/common/polylib.c | 361 ++++++++++++++++++ tools/quake3/common/polylib.h | 17 + tools/quake3/q3map2/brush.c | 140 +++++++ tools/quake3/q3map2/map.c | 212 ++++++++-- tools/quake3/q3map2/q3map2.h | 29 ++ 9 files changed, 1065 insertions(+), 59 deletions(-) create mode 100644 regression_tests/q3map2/disappearing_sliver3/NOTES.txt diff --git a/libs/mathlib.h b/libs/mathlib.h index d206d330..62edd115 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 #ifdef __cplusplus @@ -40,6 +41,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 @@ -83,6 +90,9 @@ extern const vec3_t g_vec3_axis_z; qboolean VectorCompare (const vec3_t v1, const vec3_t v2); +qboolean VectorIsOnAxis(vec3_t v); +qboolean VectorIsOnAxialPlane(vec3_t v); + vec_t VectorLength(const vec3_t v); void VectorMA( const vec3_t va, vec_t scale, const vec3_t vb, vec3_t vc ); @@ -419,6 +429,50 @@ vec_t ray_intersect_plane(const ray_t* ray, const vec3_t normal, vec_t dist); int plane_intersect_planes(const vec4_t plane1, const vec4_t plane2, const vec4_t plane3, vec3_t intersection); + +//////////////////////////////////////////////////////////////////////////////// +// 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 c549cb10..1286ec3a 100644 --- a/libs/mathlib/mathlib.c +++ b/libs/mathlib/mathlib.c @@ -30,6 +30,54 @@ const vec3_t g_vec3_axis_x = { 1, 0, 0, }; const vec3_t g_vec3_axis_y = { 0, 1, 0, }; const vec3_t g_vec3_axis_z = { 0, 0, 1, }; +/* +================ +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 @@ -119,21 +167,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 ) { @@ -584,3 +641,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/regression_tests/q3map2/disappearing_sliver/winding_logging.patch b/regression_tests/q3map2/disappearing_sliver/winding_logging.patch index c697eb92..a1babe91 100644 --- a/regression_tests/q3map2/disappearing_sliver/winding_logging.patch +++ b/regression_tests/q3map2/disappearing_sliver/winding_logging.patch @@ -1,18 +1,17 @@ Index: tools/quake3/q3map2/brush.c =================================================================== ---- tools/quake3/q3map2/brush.c (revision 371) +--- tools/quake3/q3map2/brush.c (revision 391) +++ tools/quake3/q3map2/brush.c (working copy) -@@ -356,17 +356,29 @@ - winding_t *w; +@@ -421,10 +421,16 @@ side_t *side; plane_t *plane; -+ + +- + static int brushord = -1; + brushord++; + + Sys_Printf("In CreateBrushWindings() for brush %i\n", brushord); - -- ++ /* walk the list of brush sides */ for( i = 0; i < brush->numsides; i++ ) { @@ -21,33 +20,46 @@ Index: tools/quake3/q3map2/brush.c /* get side and plane */ side = &brush->sides[ i ]; plane = &mapplanes[ side->planenum ]; - - /* make huge winding */ +@@ -435,7 +441,13 @@ + #else w = BaseWindingForPlane( plane->normal, plane->dist ); + #endif +- + + Sys_Printf(" Before clipping we have:\n"); + int z; + for (z = 0; z < w->numpoints; z++) { + Sys_Printf(" (%.8f %.8f %.8f)\n", w->p[z][0], w->p[z][1], w->p[z][2]); + } - ++ /* walk the list of brush sides */ for( j = 0; j < brush->numsides && w != NULL; j++ ) -@@ -379,6 +391,11 @@ - continue; - plane = &mapplanes[ brush->sides[ j ].planenum ^ 1 ]; + { +@@ -451,7 +463,20 @@ + #else ChopWindingInPlace( &w, plane->normal, plane->dist, 0 ); // CLIP_EPSILON ); + #endif +- + + Sys_Printf(" After clipping w/ side %i we have:\n", j); -+ for (z = 0; z < w->numpoints; z++) { -+ Sys_Printf(" (%.8f %.8f %.8f)\n", w->p[z][0], w->p[z][1], w->p[z][2]); ++ if (w) ++ { ++ for (z = 0; z < w->numpoints; z++) ++ { ++ Sys_Printf(" (%.8f %.8f %.8f)\n", w->p[z][0], w->p[z][1], w->p[z][2]); ++ } + } - ++ else ++ { ++ Sys_Printf(" winding is NULL\n"); ++ } ++ /* ydnar: fix broken windings that would generate trifans */ - FixWinding( w ); + #if EXPERIMENTAL_HIGH_PRECISION_MATH_Q3MAP2_FIXES + FixWindingAccu(w); Index: tools/quake3/q3map2/map.c =================================================================== ---- tools/quake3/q3map2/map.c (revision 371) +--- tools/quake3/q3map2/map.c (revision 391) +++ tools/quake3/q3map2/map.c (working copy) @@ -803,7 +803,11 @@ char shader[ MAX_QPATH ]; @@ -74,7 +86,7 @@ Index: tools/quake3/q3map2/map.c if( !GetToken( qtrue ) ) break; if( !strcmp( token, "}" ) ) -@@ -917,6 +924,10 @@ +@@ -917,7 +924,16 @@ } /* find the plane number */ @@ -83,5 +95,11 @@ Index: tools/quake3/q3map2/map.c + Sys_Printf(" (%f %f %f)\n", planePoints[1][0], planePoints[1][1], planePoints[1][2]); + Sys_Printf(" (%f %f %f)\n", planePoints[2][0], planePoints[2][1], planePoints[2][2]); planenum = MapPlaneFromPoints( planePoints ); ++ Sys_Printf(" normal: (%.10f %.10f %.10f)\n", ++ mapplanes[planenum].normal[0], ++ mapplanes[planenum].normal[1], ++ mapplanes[planenum].normal[2]); ++ Sys_Printf(" dist: %.10f\n", mapplanes[planenum].dist); side->planenum = planenum; + /* bp: get the texture mapping for this texturedef / plane combination */ diff --git a/regression_tests/q3map2/disappearing_sliver3/NOTES.txt b/regression_tests/q3map2/disappearing_sliver3/NOTES.txt new file mode 100644 index 00000000..aebd534f --- /dev/null +++ b/regression_tests/q3map2/disappearing_sliver3/NOTES.txt @@ -0,0 +1,36 @@ +Random notes for Rambetter, don't expect to understand this: +============================================================ + +Brush 0 is the problem. + +Side 0 is the problem (under surf tri). +Side 1 is the +y 4-face. +Side 2 is the -x 4-face. +Side 3 is the -y 4-face. +side 4 is the +z tri. + +(6144, 16122) -> (6784, 16241) +x "climb" of side 1 is 6784 - 6144 = 640. +y "climb" of side 1 is 16241 - 16122 = 119. + +x/y "climb rate" of side 1 is 640 / 119 = 5.378151261. + +After clipping side 0 against side 1, we get +************ +**** (-262144, -33762.8125) -> (262144, 63722) +************ +The slope of that is (262144 + 262144) / (63722 + 33762.8125) = 5.378150571. + +(-262144, y) -> (6784, 16241) +So (6784 + 262144) / (16241 - y) = 640 / 119 +So y = 16241 - ((119 * (6784 + 262144)) / 640) = -33762.8 + +(6144, 16122) -> (262144, y) +So (262144 - 6144) / (y - 16122) = 640 / 119 +So y = 16122 + ((119 * (262144 - 6144)) / 640) = 63722 + +After clipping side 0 against side 1 should have +************ +**** (-262144, -33762.8) -> (262144, 63722) +************ + diff --git a/tools/quake3/common/polylib.c b/tools/quake3/common/polylib.c index fc8805e7..34dd173e 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) size_t size; winding_t *c; + if (!w) Error("CopyWinding: winding is NULL"); + c = AllocWinding (w->numpoints); size = (size_t)((winding_t *)NULL)->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 d78c780d..1f4495d9 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 d61f51be..fcb6ac7d 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 */ @@ -152,28 +153,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 } @@ -192,18 +256,70 @@ void SnapPlane( vec3_t normal, vec_t *dist, vec3_t center ) SnapPlane reenabled by namespace because of multiple reports of q3map2-crashes which were triggered by this patch. */ - // div0: ensure the point "center" stays on the plane (actually, this - // rotates the plane around the point center). - // if center lies on the plane, it is guaranteed to stay on the plane by - // this fix. - vec_t centerDist = DotProduct(normal, center); SnapNormal( normal ); - *dist += (DotProduct(normal, center) - centerDist); + + // 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; + } + } +} + /* @@ -221,16 +337,15 @@ int FindFloatPlane( vec3_t innormal, vec_t dist, int numPoints, vec3_t *points ) int pidx; plane_t *p; vec_t d; - vec3_t centerofweight; vec3_t normal; - VectorClear(centerofweight); - for(i = 0; i < numPoints; ++i) - VectorMA(centerofweight, 1.0 / numPoints, points[i], centerofweight); - - /* hash the plane */ VectorCopy(innormal, normal); - SnapPlane( normal, &dist, centerofweight ); +#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 */ @@ -251,9 +366,15 @@ int FindFloatPlane( vec3_t innormal, 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 ], p->normal ) - p->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 */ @@ -273,15 +394,12 @@ int FindFloatPlane( vec3_t innormal, vec_t dist, int numPoints, vec3_t *points ) plane_t *p; vec3_t normal; - - vec3_t centerofweight; - - VectorClear(centerofweight); - for(i = 0; i < numPoints; ++i) - VectorMA(centerofweight, 1.0 / numPoints, points[i], centerofweight); - VectorCopy(innormal, normal); - SnapPlane( normal, &dist, centerofweight ); +#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 ) ) @@ -301,6 +419,9 @@ int FindFloatPlane( vec3_t innormal, vec_t dist, int numPoints, vec3_t *points ) /* found a matching plane */ if( j >= numPoints ) 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 ); @@ -317,6 +438,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; @@ -332,6 +475,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 f01d353e..bf30def5 100644 --- a/tools/quake3/q3map2/q3map2.h +++ b/tools/quake3/q3map2/q3map2.h @@ -122,6 +122,12 @@ constants ------------------------------------------------------------------------------- */ +/* temporary hacks and tests (please keep off in SVN to prevent anyone's legacy map from screwing up) */ +/* 2011-01-10 TTimo says we should turn these on in SVN, so turning on now */ +#define EXPERIMENTAL_HIGH_PRECISION_MATH_Q3MAP2_FIXES 1 +#define EXPERIMENTAL_SNAP_NORMAL_FIX 1 +#define EXPERIMENTAL_SNAP_PLANE_FIX 1 + /* general */ #define MAX_QPATH 64 @@ -2026,8 +2032,31 @@ Q_EXTERN qboolean debugPortals Q_ASSIGN( qfalse ); Q_EXTERN qboolean lightmapTriangleCheck Q_ASSIGN(qfalse); Q_EXTERN qboolean lightmapExtraVisClusterNudge Q_ASSIGN(qfalse); Q_EXTERN qboolean lightmapFill 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 */ -- 2.39.2