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