tools/quake3/common/polylib.c

   1 /*
   2    Copyright (C) 1999-2007 id Software, Inc. and contributors.
   3    For a list of contributors, see the accompanying CONTRIBUTORS file.
   4
   5    This file is part of GtkRadiant.
   6
   7    GtkRadiant is free software; you can redistribute it and/or modify
   8    it under the terms of the GNU General Public License as published by
   9    the Free Software Foundation; either version 2 of the License, or
  10    (at your option) any later version.
  11
  12    GtkRadiant is distributed in the hope that it will be useful,
  13    but WITHOUT ANY WARRANTY; without even the implied warranty of
  14    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15    GNU General Public License for more details.
  16
  17    You should have received a copy of the GNU General Public License
  18    along with GtkRadiant; if not, write to the Free Software
  19    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  20  */
  21
  22 #include "cmdlib.h"
  23 #include "mathlib.h"
  24 #include "inout.h"
  25 #include "polylib.h"
  26 #include "qfiles.h"
  27
  28
  29 extern int numthreads;
  30
  31 // counters are only bumped when running single threaded,
  32 // because they are an awefull coherence problem
  33 int c_active_windings;
  34 int c_peak_windings;
  35 int c_winding_allocs;
  36 int c_winding_points;
  37
  38 #define BOGUS_RANGE WORLD_SIZE
  39
  40 void pw( winding_t *w ){
  41         int i;
  42         for ( i = 0 ; i < w->numpoints ; i++ )
  43                 Sys_Printf( "(%5.1f, %5.1f, %5.1f)\n",w->p[i][0], w->p[i][1],w->p[i][2] );
  44 }
  45
  46
  47 /*
  48    =============
  49    AllocWinding
  50    =============
  51  */
  52 winding_t   *AllocWinding( int points ){
  53         winding_t   *w;
  54         int s;
  55
  56         if ( points >= MAX_POINTS_ON_WINDING ) {
  57                 Error( "AllocWinding failed: MAX_POINTS_ON_WINDING exceeded" );
  58         }
  59
  60         if ( numthreads == 1 ) {
  61                 c_winding_allocs++;
  62                 c_winding_points += points;
  63                 c_active_windings++;
  64                 if ( c_active_windings > c_peak_windings ) {
  65                         c_peak_windings = c_active_windings;
  66                 }
  67         }
  68         s = sizeof( *w ) + ( points ? sizeof( w->p[0] ) * ( points - 1 ) : 0 );
  69         w = safe_malloc0( s );
  70         return w;
  71 }
  72
  73 /*
  74    =============
  75    AllocWindingAccu
  76    =============
  77  */
  78 winding_accu_t *AllocWindingAccu( int points ){
  79         winding_accu_t  *w;
  80         int s;
  81
  82         if ( points >= MAX_POINTS_ON_WINDING ) {
  83                 Error( "AllocWindingAccu failed: MAX_POINTS_ON_WINDING exceeded" );
  84         }
  85
  86         if ( numthreads == 1 ) {
  87                 // At the time of this writing, these statistics were not used in any way.
  88                 c_winding_allocs++;
  89                 c_winding_points += points;
  90                 c_active_windings++;
  91                 if ( c_active_windings > c_peak_windings ) {
  92                         c_peak_windings = c_active_windings;
  93                 }
  94         }
  95         s = sizeof( *w ) + ( points ? sizeof( w->p[0] ) * ( points - 1 ) : 0 );
  96         w = safe_malloc0( s );
  97         return w;
  98 }
  99
 100 /*
 101    =============
 102    FreeWinding
 103    =============
 104  */
 105 void FreeWinding( winding_t *w ){
 106         if ( !w ) {
 107                 Error( "FreeWinding: winding is NULL" );
 108         }
 109
 110         if ( *(unsigned *)w == 0xdeaddead ) {
 111                 Error( "FreeWinding: freed a freed winding" );
 112         }
 113         *(unsigned *)w = 0xdeaddead;
 114
 115         if ( numthreads == 1 ) {
 116                 c_active_windings--;
 117         }
 118         free( w );
 119 }
 120
 121 /*
 122    =============
 123    FreeWindingAccu
 124    =============
 125  */
 126 void FreeWindingAccu( winding_accu_t *w ){
 127         if ( !w ) {
 128                 Error( "FreeWindingAccu: winding is NULL" );
 129         }
 130
 131         if ( *( (unsigned *) w ) == 0xdeaddead ) {
 132                 Error( "FreeWindingAccu: freed a freed winding" );
 133         }
 134         *( (unsigned *) w ) = 0xdeaddead;
 135
 136         if ( numthreads == 1 ) {
 137                 c_active_windings--;
 138         }
 139         free( w );
 140 }
 141
 142 /*
 143    ============
 144    RemoveColinearPoints
 145    ============
 146  */
 147 int c_removed;
 148
 149 void    RemoveColinearPoints( winding_t *w ){
 150         int i, j, k;
 151         vec3_t v1, v2;
 152         int nump;
 153         vec3_t p[MAX_POINTS_ON_WINDING];
 154
 155         nump = 0;
 156         for ( i = 0 ; i < w->numpoints ; i++ )
 157         {
 158                 j = ( i + 1 ) % w->numpoints;
 159                 k = ( i + w->numpoints - 1 ) % w->numpoints;
 160                 VectorSubtract( w->p[j], w->p[i], v1 );
 161                 VectorSubtract( w->p[i], w->p[k], v2 );
 162                 VectorNormalize( v1,v1 );
 163                 VectorNormalize( v2,v2 );
 164                 if ( DotProduct( v1, v2 ) < 0.999 ) {
 165                         VectorCopy( w->p[i], p[nump] );
 166                         nump++;
 167                 }
 168         }
 169
 170         if ( nump == w->numpoints ) {
 171                 return;
 172         }
 173
 174         if ( numthreads == 1 ) {
 175                 c_removed += w->numpoints - nump;
 176         }
 177         w->numpoints = nump;
 178         memcpy( w->p, p, nump * sizeof( p[0] ) );
 179 }
 180
 181 /*
 182    ============
 183    WindingPlane
 184    ============
 185  */
 186 void WindingPlane( winding_t *w, vec3_t normal, vec_t *dist ){
 187         vec3_t v1, v2;
 188
 189         VectorSubtract( w->p[1], w->p[0], v1 );
 190         VectorSubtract( w->p[2], w->p[0], v2 );
 191         CrossProduct( v2, v1, normal );
 192         VectorNormalize( normal, normal );
 193         *dist = DotProduct( w->p[0], normal );
 194
 195 }
 196
 197 /*
 198    =============
 199    WindingArea
 200    =============
 201  */
 202 vec_t   WindingArea( winding_t *w ){
 203         int i;
 204         vec3_t d1, d2, cross;
 205         vec_t total;
 206
 207         total = 0;
 208         for ( i = 2 ; i < w->numpoints ; i++ )
 209         {
 210                 VectorSubtract( w->p[i - 1], w->p[0], d1 );
 211                 VectorSubtract( w->p[i], w->p[0], d2 );
 212                 CrossProduct( d1, d2, cross );
 213                 total += 0.5 * VectorLength( cross );
 214         }
 215         return total;
 216 }
 217
 218 void    WindingBounds( winding_t *w, vec3_t mins, vec3_t maxs ){
 219         vec_t v;
 220         int i,j;
 221
 222         mins[0] = mins[1] = mins[2] = 99999;
 223         maxs[0] = maxs[1] = maxs[2] = -99999;
 224
 225         for ( i = 0 ; i < w->numpoints ; i++ )
 226         {
 227                 for ( j = 0 ; j < 3 ; j++ )
 228                 {
 229                         v = w->p[i][j];
 230                         if ( v < mins[j] ) {
 231                                 mins[j] = v;
 232                         }
 233                         if ( v > maxs[j] ) {
 234                                 maxs[j] = v;
 235                         }
 236                 }
 237         }
 238 }
 239
 240 /*
 241    =============
 242    WindingCenter
 243    =============
 244  */
 245 void    WindingCenter( winding_t *w, vec3_t center ){
 246         int i;
 247         float scale;
 248
 249         VectorCopy( vec3_origin, center );
 250         for ( i = 0 ; i < w->numpoints ; i++ )
 251                 VectorAdd( w->p[i], center, center );
 252
 253         scale = 1.0 / w->numpoints;
 254         VectorScale( center, scale, center );
 255 }
 256
 257 /*
 258    =================
 259    BaseWindingForPlaneAccu
 260    =================
 261  */
 262 winding_accu_t *BaseWindingForPlaneAccu( vec3_t normal, vec_t dist ){
 263         // The goal in this function is to replicate the behavior of the original BaseWindingForPlane()
 264         // function (see below) but at the same time increasing accuracy substantially.
 265
 266         // The original code gave a preference for the vup vector to start out as (0, 0, 1), unless the
 267         // normal had a dominant Z value, in which case vup started out as (1, 0, 0).  After that, vup
 268         // was "bent" [along the plane defined by normal and vup] to become perpendicular to normal.
 269         // After that the vright vector was computed as the cross product of vup and normal.
 270
 271         // I'm constructing the winding polygon points in a fashion similar to the method used in the
 272         // original function.  Orientation is the same.  The size of the winding polygon, however, is
 273         // variable in this function (depending on the angle of normal), and is larger (by about a factor
 274         // of 2) than the winding polygon in the original function.
 275
 276         int x, i;
 277         vec_t max, v;
 278         vec3_accu_t vright, vup, org, normalAccu;
 279         winding_accu_t  *w;
 280
 281         // One of the components of normal must have a magnitiude greater than this value,
 282         // otherwise normal is not a unit vector.  This is a little bit of inexpensive
 283         // partial error checking we can do.
 284         max = 0.56; // 1 / sqrt(1^2 + 1^2 + 1^2) = 0.577350269
 285
 286         x = -1;
 287         for ( i = 0; i < 3; i++ ) {
 288                 v = (vec_t) fabs( normal[i] );
 289                 if ( v > max ) {
 290                         x = i;
 291                         max = v;
 292                 }
 293         }
 294         if ( x == -1 ) {
 295                 Error( "BaseWindingForPlaneAccu: no dominant axis found because normal is too short" );
 296         }
 297
 298         switch ( x ) {
 299         case 0:     // Fall through to next case.
 300         case 1:
 301                 vright[0] = (vec_accu_t) -normal[1];
 302                 vright[1] = (vec_accu_t) normal[0];
 303                 vright[2] = 0;
 304                 break;
 305         case 2:
 306                 vright[0] = 0;
 307                 vright[1] = (vec_accu_t) -normal[2];
 308                 vright[2] = (vec_accu_t) normal[1];
 309                 break;
 310         }
 311
 312         // vright and normal are now perpendicular; you can prove this by taking their
 313         // dot product and seeing that it's always exactly 0 (with no error).
 314
 315         // NOTE: vright is NOT a unit vector at this point.  vright will have length
 316         // not exceeding 1.0.  The minimum length that vright can achieve happens when,
 317         // for example, the Z and X components of the normal input vector are equal,
 318         // and when normal's Y component is zero.  In that case Z and X of the normal
 319         // vector are both approximately 0.70711.  The resulting vright vector in this
 320         // case will have a length of 0.70711.
 321
 322         // We're relying on the fact that MAX_WORLD_COORD is a power of 2 to keep
 323         // our calculation precise and relatively free of floating point error.
 324         // [However, the code will still work fine if that's not the case.]
 325         VectorScaleAccu( vright, ( (vec_accu_t) MAX_WORLD_COORD ) * 4.0, vright );
 326
 327         // At time time of this writing, MAX_WORLD_COORD was 65536 (2^16).  Therefore
 328         // the length of vright at this point is at least 185364.  In comparison, a
 329         // corner of the world at location (65536, 65536, 65536) is distance 113512
 330         // away from the origin.
 331
 332         VectorCopyRegularToAccu( normal, normalAccu );
 333         CrossProductAccu( normalAccu, vright, vup );
 334
 335         // vup now has length equal to that of vright.
 336
 337         VectorScaleAccu( normalAccu, (vec_accu_t) dist, org );
 338
 339         // org is now a point on the plane defined by normal and dist.  Furthermore,
 340         // org, vright, and vup are pairwise perpendicular.  Now, the 4 vectors
 341         // { (+-)vright + (+-)vup } have length that is at least sqrt(185364^2 + 185364^2),
 342         // which is about 262144.  That length lies outside the world, since the furthest
 343         // point in the world has distance 113512 from the origin as mentioned above.
 344         // Also, these 4 vectors are perpendicular to the org vector.  So adding them
 345         // to org will only increase their length.  Therefore the 4 points defined below
 346         // all lie outside of the world.  Furthermore, it can be easily seen that the
 347         // edges connecting these 4 points (in the winding_accu_t below) lie completely
 348         // outside the world.  sqrt(262144^2 + 262144^2)/2 = 185363, which is greater than
 349         // 113512.
 350
 351         w = AllocWindingAccu( 4 );
 352
 353         VectorSubtractAccu( org, vright, w->p[0] );
 354         VectorAddAccu( w->p[0], vup, w->p[0] );
 355
 356         VectorAddAccu( org, vright, w->p[1] );
 357         VectorAddAccu( w->p[1], vup, w->p[1] );
 358
 359         VectorAddAccu( org, vright, w->p[2] );
 360         VectorSubtractAccu( w->p[2], vup, w->p[2] );
 361
 362         VectorSubtractAccu( org, vright, w->p[3] );
 363         VectorSubtractAccu( w->p[3], vup, w->p[3] );
 364
 365         w->numpoints = 4;
 366
 367         return w;
 368 }
 369
 370 /*
 371    =================
 372    BaseWindingForPlane
 373
 374    Original BaseWindingForPlane() function that has serious accuracy problems.  Here is why.
 375    The base winding is computed as a rectangle with very large coordinates.  These coordinates
 376    are in the range 2^17 or 2^18.  "Epsilon" (meaning the distance between two adjacent numbers)
 377    at these magnitudes in 32 bit floating point world is about 0.02.  So for example, if things
 378    go badly (by bad luck), then the whole plane could be shifted by 0.02 units (its distance could
 379    be off by that much).  Then if we were to compute the winding of this plane and another of
 380    the brush's planes met this winding at a very acute angle, that error could multiply to around
 381    0.1 or more when computing the final vertex coordinates of the winding.  0.1 is a very large
 382    error, and can lead to all sorts of disappearing triangle problems.
 383    =================
 384  */
 385 winding_t *BaseWindingForPlane( vec3_t normal, vec_t dist ){
 386         int i, x;
 387         vec_t max, v;
 388         vec3_t org, vright, vup;
 389         winding_t   *w;
 390
 391 // find the major axis
 392
 393         max = -BOGUS_RANGE;
 394         x = -1;
 395         for ( i = 0 ; i < 3; i++ )
 396         {
 397                 v = fabs( normal[i] );
 398                 if ( v > max ) {
 399                         x = i;
 400                         max = v;
 401                 }
 402         }
 403         if ( x == -1 ) {
 404                 Error( "BaseWindingForPlane: no axis found" );
 405         }
 406
 407         VectorCopy( vec3_origin, vup );
 408         switch ( x )
 409         {
 410         case 0:
 411         case 1:
 412                 vup[2] = 1;
 413                 break;
 414         case 2:
 415                 vup[0] = 1;
 416                 break;
 417         }
 418
 419         v = DotProduct( vup, normal );
 420         VectorMA( vup, -v, normal, vup );
 421         VectorNormalize( vup, vup );
 422
 423         VectorScale( normal, dist, org );
 424
 425         CrossProduct( vup, normal, vright );
 426
 427         // LordHavoc: this has to use *2 because otherwise some created points may
 428         // be inside the world (think of a diagonal case), and any brush with such
 429         // points should be removed, failure to detect such cases is disasterous
 430         VectorScale( vup, MAX_WORLD_COORD * 2, vup );
 431         VectorScale( vright, MAX_WORLD_COORD * 2, vright );
 432
 433         // project a really big axis aligned box onto the plane
 434         w = AllocWinding( 4 );
 435
 436         VectorSubtract( org, vright, w->p[0] );
 437         VectorAdd( w->p[0], vup, w->p[0] );
 438
 439         VectorAdd( org, vright, w->p[1] );
 440         VectorAdd( w->p[1], vup, w->p[1] );
 441
 442         VectorAdd( org, vright, w->p[2] );
 443         VectorSubtract( w->p[2], vup, w->p[2] );
 444
 445         VectorSubtract( org, vright, w->p[3] );
 446         VectorSubtract( w->p[3], vup, w->p[3] );
 447
 448         w->numpoints = 4;
 449
 450         return w;
 451 }
 452
 453 /*
 454    ==================
 455    CopyWinding
 456    ==================
 457  */
 458 winding_t   *CopyWinding( winding_t *w ){
 459         size_t size;
 460         winding_t   *c;
 461
 462         if ( !w ) {
 463                 Error( "CopyWinding: winding is NULL" );
 464         }
 465
 466         c = AllocWinding( w->numpoints );
 467         size = (size_t)( (winding_t *)NULL )->p[w->numpoints];
 468         memcpy( c, w, size );
 469         return c;
 470 }
 471
 472 /*
 473    ==================
 474    CopyWindingAccuIncreaseSizeAndFreeOld
 475    ==================
 476  */
 477 winding_accu_t *CopyWindingAccuIncreaseSizeAndFreeOld( winding_accu_t *w ){
 478         int i;
 479         winding_accu_t  *c;
 480
 481         if ( !w ) {
 482                 Error( "CopyWindingAccuIncreaseSizeAndFreeOld: winding is NULL" );
 483         }
 484
 485         c = AllocWindingAccu( w->numpoints + 1 );
 486         c->numpoints = w->numpoints;
 487         for ( i = 0; i < c->numpoints; i++ )
 488         {
 489                 VectorCopyAccu( w->p[i], c->p[i] );
 490         }
 491         FreeWindingAccu( w );
 492         return c;
 493 }
 494
 495 /*
 496    ==================
 497    CopyWindingAccuToRegular
 498    ==================
 499  */
 500 winding_t   *CopyWindingAccuToRegular( winding_accu_t *w ){
 501         int i;
 502         winding_t   *c;
 503
 504         if ( !w ) {
 505                 Error( "CopyWindingAccuToRegular: winding is NULL" );
 506         }
 507
 508         c = AllocWinding( w->numpoints );
 509         c->numpoints = w->numpoints;
 510         for ( i = 0; i < c->numpoints; i++ )
 511         {
 512                 VectorCopyAccuToRegular( w->p[i], c->p[i] );
 513         }
 514         return c;
 515 }
 516
 517 /*
 518    ==================
 519    ReverseWinding
 520    ==================
 521  */
 522 winding_t   *ReverseWinding( winding_t *w ){
 523         int i;
 524         winding_t   *c;
 525
 526         c = AllocWinding( w->numpoints );
 527         for ( i = 0 ; i < w->numpoints ; i++ )
 528         {
 529                 VectorCopy( w->p[w->numpoints - 1 - i], c->p[i] );
 530         }
 531         c->numpoints = w->numpoints;
 532         return c;
 533 }
 534
 535
 536 /*
 537    =============
 538    ClipWindingEpsilon
 539    =============
 540  */
 541 void    ClipWindingEpsilonStrict( winding_t *in, vec3_t normal, vec_t dist,
 542                                                         vec_t epsilon, winding_t **front, winding_t **back ){
 543         vec_t dists[MAX_POINTS_ON_WINDING + 4];
 544         int sides[MAX_POINTS_ON_WINDING + 4];
 545         int counts[3];
 546         static vec_t dot;           // VC 4.2 optimizer bug if not static
 547         int i, j;
 548         vec_t   *p1, *p2;
 549         vec3_t mid;
 550         winding_t   *f, *b;
 551         int maxpts;
 552
 553         counts[0] = counts[1] = counts[2] = 0;
 554
 555 // determine sides for each point
 556         for ( i = 0 ; i < in->numpoints ; i++ )
 557         {
 558
 559                 dot = DotProduct( in->p[i], normal );
 560                 dot -= dist;
 561                 dists[i] = dot;
 562                 if ( dot > epsilon ) {
 563                         sides[i] = SIDE_FRONT;
 564                 }
 565                 else if ( dot < -epsilon ) {
 566                         sides[i] = SIDE_BACK;
 567                 }
 568                 else
 569                 {
 570                         sides[i] = SIDE_ON;
 571                 }
 572                 counts[sides[i]]++;
 573         }
 574         sides[i] = sides[0];
 575         dists[i] = dists[0];
 576
 577         *front = *back = NULL;
 578
 579         if ( !counts[0] && !counts[1] ) {
 580                 return;
 581         }
 582         if ( !counts[0] ) {
 583                 *back = CopyWinding( in );
 584                 return;
 585         }
 586         if ( !counts[1] ) {
 587                 *front = CopyWinding( in );
 588                 return;
 589         }
 590
 591         maxpts = in->numpoints + 4;   // cant use counts[0]+2 because
 592                                       // of fp grouping errors
 593
 594         *front = f = AllocWinding( maxpts );
 595         *back = b = AllocWinding( maxpts );
 596
 597         for ( i = 0 ; i < in->numpoints ; i++ )
 598         {
 599                 p1 = in->p[i];
 600
 601                 if ( sides[i] == SIDE_ON ) {
 602                         VectorCopy( p1, f->p[f->numpoints] );
 603                         f->numpoints++;
 604                         VectorCopy( p1, b->p[b->numpoints] );
 605                         b->numpoints++;
 606                         continue;
 607                 }
 608
 609                 if ( sides[i] == SIDE_FRONT ) {
 610                         VectorCopy( p1, f->p[f->numpoints] );
 611                         f->numpoints++;
 612                 }
 613                 if ( sides[i] == SIDE_BACK ) {
 614                         VectorCopy( p1, b->p[b->numpoints] );
 615                         b->numpoints++;
 616                 }
 617
 618                 if ( sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i] ) {
 619                         continue;
 620                 }
 621
 622                 // generate a split point
 623                 p2 = in->p[( i + 1 ) % in->numpoints];
 624
 625                 dot = dists[i] / ( dists[i] - dists[i + 1] );
 626                 for ( j = 0 ; j < 3 ; j++ )
 627                 {   // avoid round off error when possible
 628                         if ( normal[j] == 1 ) {
 629                                 mid[j] = dist;
 630                         }
 631                         else if ( normal[j] == -1 ) {
 632                                 mid[j] = -dist;
 633                         }
 634                         else{
 635                                 mid[j] = p1[j] + dot * ( p2[j] - p1[j] );
 636                         }
 637                 }
 638
 639                 VectorCopy( mid, f->p[f->numpoints] );
 640                 f->numpoints++;
 641                 VectorCopy( mid, b->p[b->numpoints] );
 642                 b->numpoints++;
 643         }
 644
 645         if ( f->numpoints > maxpts || b->numpoints > maxpts ) {
 646                 Error( "ClipWinding: points exceeded estimate" );
 647         }
 648         if ( f->numpoints > MAX_POINTS_ON_WINDING || b->numpoints > MAX_POINTS_ON_WINDING ) {
 649                 Error( "ClipWinding: MAX_POINTS_ON_WINDING" );
 650         }
 651 }
 652
 653 void    ClipWindingEpsilon( winding_t *in, vec3_t normal, vec_t dist,
 654                                                         vec_t epsilon, winding_t **front, winding_t **back ){
 655         ClipWindingEpsilonStrict( in, normal, dist, epsilon, front, back );
 656         /* apparently most code expects that in the winding-on-plane case, the back winding is the original winding */
 657         if ( !*front && !*back ) {
 658                 *back = CopyWinding( in );
 659         }
 660 }
 661
 662
 663 /*
 664    =============
 665    ChopWindingInPlaceAccu
 666    =============
 667  */
 668 void ChopWindingInPlaceAccu( winding_accu_t **inout, vec3_t normal, vec_t dist, vec_t crudeEpsilon ){
 669         vec_accu_t fineEpsilon;
 670         winding_accu_t  *in;
 671         int counts[3];
 672         int i, j;
 673         vec_accu_t dists[MAX_POINTS_ON_WINDING + 1];
 674         int sides[MAX_POINTS_ON_WINDING + 1];
 675         int maxpts;
 676         winding_accu_t  *f;
 677         vec_accu_t  *p1, *p2;
 678         vec_accu_t w;
 679         vec3_accu_t mid, normalAccu;
 680
 681         // We require at least a very small epsilon.  It's a good idea for several reasons.
 682         // First, we will be dividing by a potentially very small distance below.  We don't
 683         // want that distance to be too small; otherwise, things "blow up" with little accuracy
 684         // due to the division.  (After a second look, the value w below is in range (0,1), but
 685         // graininess problem remains.)  Second, Having minimum epsilon also prevents the following
 686         // situation.  Say for example we have a perfect octagon defined by the input winding.
 687         // Say our chopping plane (defined by normal and dist) is essentially the same plane
 688         // that the octagon is sitting on.  Well, due to rounding errors, it may be that point
 689         // 1 of the octagon might be in front, point 2 might be in back, point 3 might be in
 690         // front, point 4 might be in back, and so on.  So we could end up with a very ugly-
 691         // looking chopped winding, and this might be undesirable, and would at least lead to
 692         // a possible exhaustion of MAX_POINTS_ON_WINDING.  It's better to assume that points
 693         // very very close to the plane are on the plane, using an infinitesimal epsilon amount.
 694
 695         // Now, the original ChopWindingInPlace() function used a vec_t-based winding_t.
 696         // So this minimum epsilon is quite similar to casting the higher resolution numbers to
 697         // the lower resolution and comparing them in the lower resolution mode.  We explicitly
 698         // choose the minimum epsilon as something around the vec_t epsilon of one because we
 699         // want the resolution of vec_accu_t to have a large resolution around the epsilon.
 700         // Some of that leftover resolution even goes away after we scale to points far away.
 701
 702         // Here is a further discussion regarding the choice of smallestEpsilonAllowed.
 703         // In the 32 float world (we can assume vec_t is that), the "epsilon around 1.0" is
 704         // 0.00000011921.  In the 64 bit float world (we can assume vec_accu_t is that), the
 705         // "epsilon around 1.0" is 0.00000000000000022204.  (By the way these two epsilons
 706         // are defined as VEC_SMALLEST_EPSILON_AROUND_ONE VEC_ACCU_SMALLEST_EPSILON_AROUND_ONE
 707         // respectively.)  If you divide the first by the second, you get approximately
 708         // 536,885,246.  Dividing that number by 200,000 (a typical base winding coordinate)
 709         // gives 2684.  So in other words, if our smallestEpsilonAllowed was chosen as exactly
 710         // VEC_SMALLEST_EPSILON_AROUND_ONE, you would be guaranteed at least 2000 "ticks" in
 711         // 64-bit land inside of the epsilon for all numbers we're dealing with.
 712
 713         static const vec_accu_t smallestEpsilonAllowed = ( (vec_accu_t) VEC_SMALLEST_EPSILON_AROUND_ONE ) * 0.5;
 714         if ( crudeEpsilon < smallestEpsilonAllowed ) {
 715                 fineEpsilon = smallestEpsilonAllowed;
 716         }
 717         else{fineEpsilon = (vec_accu_t) crudeEpsilon; }
 718
 719         in = *inout;
 720         counts[0] = counts[1] = counts[2] = 0;
 721         VectorCopyRegularToAccu( normal, normalAccu );
 722
 723         for ( i = 0; i < in->numpoints; i++ )
 724         {
 725                 dists[i] = DotProductAccu( in->p[i], normalAccu ) - dist;
 726                 if ( dists[i] > fineEpsilon ) {
 727                         sides[i] = SIDE_FRONT;
 728                 }
 729                 else if ( dists[i] < -fineEpsilon ) {
 730                         sides[i] = SIDE_BACK;
 731                 }
 732                 else{sides[i] = SIDE_ON; }
 733                 counts[sides[i]]++;
 734         }
 735         sides[i] = sides[0];
 736         dists[i] = dists[0];
 737
 738         // I'm wondering if whatever code that handles duplicate planes is robust enough
 739         // that we never get a case where two nearly equal planes result in 2 NULL windings
 740         // due to the 'if' statement below.  TODO: Investigate this.
 741         if ( !counts[SIDE_FRONT] ) {
 742                 FreeWindingAccu( in );
 743                 *inout = NULL;
 744                 return;
 745         }
 746         if ( !counts[SIDE_BACK] ) {
 747                 return; // Winding is unmodified.
 748         }
 749
 750         // NOTE: The least number of points that a winding can have at this point is 2.
 751         // In that case, one point is SIDE_FRONT and the other is SIDE_BACK.
 752
 753         maxpts = counts[SIDE_FRONT] + 2; // We dynamically expand if this is too small.
 754         f = AllocWindingAccu( maxpts );
 755
 756         for ( i = 0; i < in->numpoints; i++ )
 757         {
 758                 p1 = in->p[i];
 759
 760                 if ( sides[i] == SIDE_ON || sides[i] == SIDE_FRONT ) {
 761                         if ( f->numpoints >= MAX_POINTS_ON_WINDING ) {
 762                                 Error( "ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING" );
 763                         }
 764                         if ( f->numpoints >= maxpts ) { // This will probably never happen.
 765                                 Sys_FPrintf( SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n" );
 766                                 f = CopyWindingAccuIncreaseSizeAndFreeOld( f );
 767                                 maxpts++;
 768                         }
 769                         VectorCopyAccu( p1, f->p[f->numpoints] );
 770                         f->numpoints++;
 771                         if ( sides[i] == SIDE_ON ) {
 772                                 continue;
 773                         }
 774                 }
 775                 if ( sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i] ) {
 776                         continue;
 777                 }
 778
 779                 // Generate a split point.
 780                 p2 = in->p[( ( i + 1 ) == in->numpoints ) ? 0 : ( i + 1 )];
 781
 782                 // The divisor's absolute value is greater than the dividend's absolute value.
 783                 // w is in the range (0,1).
 784                 w = dists[i] / ( dists[i] - dists[i + 1] );
 785
 786                 for ( j = 0; j < 3; j++ )
 787                 {
 788                         // Avoid round-off error when possible.  Check axis-aligned normal.
 789                         if ( normal[j] == 1 ) {
 790                                 mid[j] = dist;
 791                         }
 792                         else if ( normal[j] == -1 ) {
 793                                 mid[j] = -dist;
 794                         }
 795                         else{mid[j] = p1[j] + ( w * ( p2[j] - p1[j] ) ); }
 796                 }
 797                 if ( f->numpoints >= MAX_POINTS_ON_WINDING ) {
 798                         Error( "ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING" );
 799                 }
 800                 if ( f->numpoints >= maxpts ) { // This will probably never happen.
 801                         Sys_FPrintf( SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n" );
 802                         f = CopyWindingAccuIncreaseSizeAndFreeOld( f );
 803                         maxpts++;
 804                 }
 805                 VectorCopyAccu( mid, f->p[f->numpoints] );
 806                 f->numpoints++;
 807         }
 808
 809         FreeWindingAccu( in );
 810         *inout = f;
 811 }
 812
 813 /*
 814    =============
 815    ChopWindingInPlace
 816    =============
 817  */
 818 void ChopWindingInPlace( winding_t **inout, vec3_t normal, vec_t dist, vec_t epsilon ){
 819         winding_t   *in;
 820         vec_t dists[MAX_POINTS_ON_WINDING + 4];
 821         int sides[MAX_POINTS_ON_WINDING + 4];
 822         int counts[3];
 823         static vec_t dot;           // VC 4.2 optimizer bug if not static
 824         int i, j;
 825         vec_t   *p1, *p2;
 826         vec3_t mid;
 827         winding_t   *f;
 828         int maxpts;
 829
 830         in = *inout;
 831         counts[0] = counts[1] = counts[2] = 0;
 832
 833 // determine sides for each point
 834         for ( i = 0 ; i < in->numpoints ; i++ )
 835         {
 836                 dot = DotProduct( in->p[i], normal );
 837                 dot -= dist;
 838                 dists[i] = dot;
 839                 if ( dot > epsilon ) {
 840                         sides[i] = SIDE_FRONT;
 841                 }
 842                 else if ( dot < -epsilon ) {
 843                         sides[i] = SIDE_BACK;
 844                 }
 845                 else
 846                 {
 847                         sides[i] = SIDE_ON;
 848                 }
 849                 counts[sides[i]]++;
 850         }
 851         sides[i] = sides[0];
 852         dists[i] = dists[0];
 853
 854         if ( !counts[0] ) {
 855                 FreeWinding( in );
 856                 *inout = NULL;
 857                 return;
 858         }
 859         if ( !counts[1] ) {
 860                 return;     // inout stays the same
 861
 862         }
 863         maxpts = in->numpoints + 4;   // cant use counts[0]+2 because
 864                                       // of fp grouping errors
 865
 866         f = AllocWinding( maxpts );
 867
 868         for ( i = 0 ; i < in->numpoints ; i++ )
 869         {
 870                 p1 = in->p[i];
 871
 872                 if ( sides[i] == SIDE_ON ) {
 873                         VectorCopy( p1, f->p[f->numpoints] );
 874                         f->numpoints++;
 875                         continue;
 876                 }
 877
 878                 if ( sides[i] == SIDE_FRONT ) {
 879                         VectorCopy( p1, f->p[f->numpoints] );
 880                         f->numpoints++;
 881                 }
 882
 883                 if ( sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i] ) {
 884                         continue;
 885                 }
 886
 887                 // generate a split point
 888                 p2 = in->p[( i + 1 ) % in->numpoints];
 889
 890                 dot = dists[i] / ( dists[i] - dists[i + 1] );
 891                 for ( j = 0 ; j < 3 ; j++ )
 892                 {   // avoid round off error when possible
 893                         if ( normal[j] == 1 ) {
 894                                 mid[j] = dist;
 895                         }
 896                         else if ( normal[j] == -1 ) {
 897                                 mid[j] = -dist;
 898                         }
 899                         else{
 900                                 mid[j] = p1[j] + dot * ( p2[j] - p1[j] );
 901                         }
 902                 }
 903
 904                 VectorCopy( mid, f->p[f->numpoints] );
 905                 f->numpoints++;
 906         }
 907
 908         if ( f->numpoints > maxpts ) {
 909                 Error( "ClipWinding: points exceeded estimate" );
 910         }
 911         if ( f->numpoints > MAX_POINTS_ON_WINDING ) {
 912                 Error( "ClipWinding: MAX_POINTS_ON_WINDING" );
 913         }
 914
 915         FreeWinding( in );
 916         *inout = f;
 917 }
 918
 919
 920 /*
 921    =================
 922    ChopWinding
 923
 924    Returns the fragment of in that is on the front side
 925    of the cliping plane.  The original is freed.
 926    =================
 927  */
 928 winding_t   *ChopWinding( winding_t *in, vec3_t normal, vec_t dist ){
 929         winding_t   *f, *b;
 930
 931         ClipWindingEpsilon( in, normal, dist, ON_EPSILON, &f, &b );
 932         FreeWinding( in );
 933         if ( b ) {
 934                 FreeWinding( b );
 935         }
 936         return f;
 937 }
 938
 939
 940 /*
 941    =================
 942    CheckWinding
 943
 944    =================
 945  */
 946 void CheckWinding( winding_t *w ){
 947         int i, j;
 948         vec_t   *p1, *p2;
 949         vec_t d, edgedist;
 950         vec3_t dir, edgenormal, facenormal;
 951         vec_t area;
 952         vec_t facedist;
 953
 954         if ( w->numpoints < 3 ) {
 955                 Error( "CheckWinding: %i points",w->numpoints );
 956         }
 957
 958         area = WindingArea( w );
 959         if ( area < 1 ) {
 960                 Error( "CheckWinding: %f area", area );
 961         }
 962
 963         WindingPlane( w, facenormal, &facedist );
 964
 965         for ( i = 0 ; i < w->numpoints ; i++ )
 966         {
 967                 p1 = w->p[i];
 968
 969                 for ( j = 0 ; j < 3 ; j++ )
 970                         if ( p1[j] > MAX_WORLD_COORD || p1[j] < MIN_WORLD_COORD ) {
 971                                 Error( "CheckFace: MAX_WORLD_COORD exceeded: %f",p1[j] );
 972                         }
 973
 974                 j = i + 1 == w->numpoints ? 0 : i + 1;
 975
 976                 // check the point is on the face plane
 977                 d = DotProduct( p1, facenormal ) - facedist;
 978                 if ( d < -ON_EPSILON || d > ON_EPSILON ) {
 979                         Error( "CheckWinding: point off plane" );
 980                 }
 981
 982                 // check the edge isnt degenerate
 983                 p2 = w->p[j];
 984                 VectorSubtract( p2, p1, dir );
 985
 986                 if ( VectorLength( dir ) < ON_EPSILON ) {
 987                         Error( "CheckWinding: degenerate edge" );
 988                 }
 989
 990                 CrossProduct( facenormal, dir, edgenormal );
 991                 VectorNormalize( edgenormal, edgenormal );
 992                 edgedist = DotProduct( p1, edgenormal );
 993                 edgedist += ON_EPSILON;
 994
 995                 // all other points must be on front side
 996                 for ( j = 0 ; j < w->numpoints ; j++ )
 997                 {
 998                         if ( j == i ) {
 999                                 continue;
1000                         }
1001                         d = DotProduct( w->p[j], edgenormal );
1002                         if ( d > edgedist ) {
1003                                 Error( "CheckWinding: non-convex" );
1004                         }
1005                 }
1006         }
1007 }
1008
1009
1010 /*
1011    ============
1012    WindingOnPlaneSide
1013    ============
1014  */
1015 int     WindingOnPlaneSide( winding_t *w, vec3_t normal, vec_t dist ){
1016         qboolean front, back;
1017         int i;
1018         vec_t d;
1019
1020         front = qfalse;
1021         back = qfalse;
1022         for ( i = 0 ; i < w->numpoints ; i++ )
1023         {
1024                 d = DotProduct( w->p[i], normal ) - dist;
1025                 if ( d < -ON_EPSILON ) {
1026                         if ( front ) {
1027                                 return SIDE_CROSS;
1028                         }
1029                         back = qtrue;
1030                         continue;
1031                 }
1032                 if ( d > ON_EPSILON ) {
1033                         if ( back ) {
1034                                 return SIDE_CROSS;
1035                         }
1036                         front = qtrue;
1037                         continue;
1038                 }
1039         }
1040
1041         if ( back ) {
1042                 return SIDE_BACK;
1043         }
1044         if ( front ) {
1045                 return SIDE_FRONT;
1046         }
1047         return SIDE_ON;
1048 }
1049
1050
1051 /*
1052    =================
1053    AddWindingToConvexHull
1054
1055    Both w and *hull are on the same plane
1056    =================
1057  */
1058 #define MAX_HULL_POINTS     128
1059 void    AddWindingToConvexHull( winding_t *w, winding_t **hull, vec3_t normal ) {
1060         int i, j, k;
1061         float       *p, *copy;
1062         vec3_t dir;
1063         float d;
1064         int numHullPoints, numNew;
1065         vec3_t hullPoints[MAX_HULL_POINTS];
1066         vec3_t newHullPoints[MAX_HULL_POINTS];
1067         vec3_t hullDirs[MAX_HULL_POINTS];
1068         qboolean hullSide[MAX_HULL_POINTS];
1069         qboolean outside;
1070
1071         if ( !*hull ) {
1072                 *hull = CopyWinding( w );
1073                 return;
1074         }
1075
1076         numHullPoints = ( *hull )->numpoints;
1077         memcpy( hullPoints, ( *hull )->p, numHullPoints * sizeof( vec3_t ) );
1078
1079         for ( i = 0 ; i < w->numpoints ; i++ ) {
1080                 p = w->p[i];
1081
1082                 // calculate hull side vectors
1083                 for ( j = 0 ; j < numHullPoints ; j++ ) {
1084                         k = ( j + 1 ) % numHullPoints;
1085
1086                         VectorSubtract( hullPoints[k], hullPoints[j], dir );
1087                         VectorNormalize( dir, dir );
1088                         CrossProduct( normal, dir, hullDirs[j] );
1089                 }
1090
1091                 outside = qfalse;
1092                 for ( j = 0 ; j < numHullPoints ; j++ ) {
1093                         VectorSubtract( p, hullPoints[j], dir );
1094                         d = DotProduct( dir, hullDirs[j] );
1095                         if ( d >= ON_EPSILON ) {
1096                                 outside = qtrue;
1097                         }
1098                         if ( d >= -ON_EPSILON ) {
1099                                 hullSide[j] = qtrue;
1100                         }
1101                         else {
1102                                 hullSide[j] = qfalse;
1103                         }
1104                 }
1105
1106                 // if the point is effectively inside, do nothing
1107                 if ( !outside ) {
1108                         continue;
1109                 }
1110
1111                 // find the back side to front side transition
1112                 for ( j = 0 ; j < numHullPoints ; j++ ) {
1113                         if ( !hullSide[ j % numHullPoints ] && hullSide[ ( j + 1 ) % numHullPoints ] ) {
1114                                 break;
1115                         }
1116                 }
1117                 if ( j == numHullPoints ) {
1118                         continue;
1119                 }
1120
1121                 // insert the point here
1122                 VectorCopy( p, newHullPoints[0] );
1123                 numNew = 1;
1124
1125                 // copy over all points that aren't double fronts
1126                 j = ( j + 1 ) % numHullPoints;
1127                 for ( k = 0 ; k < numHullPoints ; k++ ) {
1128                         if ( hullSide[ ( j + k ) % numHullPoints ] && hullSide[ ( j + k + 1 ) % numHullPoints ] ) {
1129                                 continue;
1130                         }
1131                         copy = hullPoints[ ( j + k + 1 ) % numHullPoints ];
1132                         VectorCopy( copy, newHullPoints[numNew] );
1133                         numNew++;
1134                 }
1135
1136                 numHullPoints = numNew;
1137                 memcpy( hullPoints, newHullPoints, numHullPoints * sizeof( vec3_t ) );
1138         }
1139
1140         FreeWinding( *hull );
1141         w = AllocWinding( numHullPoints );
1142         w->numpoints = numHullPoints;
1143         *hull = w;
1144         memcpy( w->p, hullPoints, numHullPoints * sizeof( vec3_t ) );
1145 }