::zerowing-base=422
authorRudolf Polzer <divVerent@xonotic.org>
Tue, 11 Jan 2011 13:39:55 +0000 (14:39 +0100)
committerRudolf Polzer <divVerent@xonotic.org>
Tue, 11 Jan 2011 13:39:55 +0000 (14:39 +0100)
libs/mathlib.h
libs/mathlib/mathlib.c
regression_tests/q3map2/disappearing_sliver/winding_logging.patch
regression_tests/q3map2/disappearing_sliver3/NOTES.txt [new file with mode: 0644]
tools/quake3/common/polylib.c
tools/quake3/common/polylib.h
tools/quake3/q3map2/brush.c
tools/quake3/q3map2/map.c
tools/quake3/q3map2/q3map2.h

index d206d33..62edd11 100644 (file)
@@ -24,6 +24,7 @@ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 
 // mathlib.h
 #include <math.h>
+#include <float.h>
 
 #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
index c549cb1..1286ec3 100644 (file)
@@ -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;
+}
index c697eb9..a1babe9 100644 (file)
@@ -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 (file)
index 0000000..aebd534
--- /dev/null
@@ -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)
+************
+
index fc8805e..34dd173 100644 (file)
@@ -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
index d78c780..1f4495d 100644 (file)
@@ -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);
index 2517ca6..1275408 100644 (file)
@@ -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 ; 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]);
index d61f51b..fcb6ac7 100644 (file)
@@ -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
 }
 
 
index f01d353..bf30def 100644 (file)
@@ -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 */