tools/quake3/common/polylib.c

   1 /*
   2    Copyright (C) 1999-2006 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( *w ) + ( points ? sizeof( w->p[0] ) * ( points - 1 ) : 0 );
  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( *w ) + ( points ? sizeof( w->p[0] ) * ( points - 1 ) : 0 );
  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         size_t size;
 463         winding_t   *c;
 464
 465         if ( !w ) {
 466                 Error( "CopyWinding: winding is NULL" );
 467         }
 468
 469         c = AllocWinding( w->numpoints );
 470         size = (size_t)( (winding_t *)NULL )->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    ClipWindingEpsilonStrict( 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] && !counts[1] ) {
 583                 return;
 584         }
 585         if ( !counts[0] ) {
 586                 *back = CopyWinding( in );
 587                 return;
 588         }
 589         if ( !counts[1] ) {
 590                 *front = CopyWinding( in );
 591                 return;
 592         }
 593
 594         maxpts = in->numpoints + 4;   // cant use counts[0]+2 because
 595                                       // of fp grouping errors
 596
 597         *front = f = AllocWinding( maxpts );
 598         *back = b = AllocWinding( maxpts );
 599
 600         for ( i = 0 ; i < in->numpoints ; i++ )
 601         {
 602                 p1 = in->p[i];
 603
 604                 if ( sides[i] == SIDE_ON ) {
 605                         VectorCopy( p1, f->p[f->numpoints] );
 606                         f->numpoints++;
 607                         VectorCopy( p1, b->p[b->numpoints] );
 608                         b->numpoints++;
 609                         continue;
 610                 }
 611
 612                 if ( sides[i] == SIDE_FRONT ) {
 613                         VectorCopy( p1, f->p[f->numpoints] );
 614                         f->numpoints++;
 615                 }
 616                 if ( sides[i] == SIDE_BACK ) {
 617                         VectorCopy( p1, b->p[b->numpoints] );
 618                         b->numpoints++;
 619                 }
 620
 621                 if ( sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i] ) {
 622                         continue;
 623                 }
 624
 625                 // generate a split point
 626                 p2 = in->p[( i + 1 ) % in->numpoints];
 627
 628                 dot = dists[i] / ( dists[i] - dists[i + 1] );
 629                 for ( j = 0 ; j < 3 ; j++ )
 630                 {   // avoid round off error when possible
 631                         if ( normal[j] == 1 ) {
 632                                 mid[j] = dist;
 633                         }
 634                         else if ( normal[j] == -1 ) {
 635                                 mid[j] = -dist;
 636                         }
 637                         else{
 638                                 mid[j] = p1[j] + dot * ( p2[j] - p1[j] );
 639                         }
 640                 }
 641
 642                 VectorCopy( mid, f->p[f->numpoints] );
 643                 f->numpoints++;
 644                 VectorCopy( mid, b->p[b->numpoints] );
 645                 b->numpoints++;
 646         }
 647
 648         if ( f->numpoints > maxpts || b->numpoints > maxpts ) {
 649                 Error( "ClipWinding: points exceeded estimate" );
 650         }
 651         if ( f->numpoints > MAX_POINTS_ON_WINDING || b->numpoints > MAX_POINTS_ON_WINDING ) {
 652                 Error( "ClipWinding: MAX_POINTS_ON_WINDING" );
 653         }
 654 }
 655
 656 void    ClipWindingEpsilon( winding_t *in, vec3_t normal, vec_t dist,
 657                                                         vec_t epsilon, winding_t **front, winding_t **back ){
 658         ClipWindingEpsilonStrict( in, normal, dist, epsilon, front, back );
 659         /* apparently most code expects that in the winding-on-plane case, the back winding is the original winding */
 660         if ( !*front && !*back ) {
 661                 *back = CopyWinding( in );
 662         }
 663 }
 664
 665
 666 /*
 667    =============
 668    ChopWindingInPlaceAccu
 669    =============
 670  */
 671 void ChopWindingInPlaceAccu( winding_accu_t **inout, vec3_t normal, vec_t dist, vec_t crudeEpsilon ){
 672         vec_accu_t fineEpsilon;
 673         winding_accu_t  *in;
 674         int counts[3];
 675         int i, j;
 676         vec_accu_t dists[MAX_POINTS_ON_WINDING + 1];
 677         int sides[MAX_POINTS_ON_WINDING + 1];
 678         int maxpts;
 679         winding_accu_t  *f;
 680         vec_accu_t  *p1, *p2;
 681         vec_accu_t w;
 682         vec3_accu_t mid, normalAccu;
 683
 684         // We require at least a very small epsilon.  It's a good idea for several reasons.
 685         // First, we will be dividing by a potentially very small distance below.  We don't
 686         // want that distance to be too small; otherwise, things "blow up" with little accuracy
 687         // due to the division.  (After a second look, the value w below is in range (0,1), but
 688         // graininess problem remains.)  Second, Having minimum epsilon also prevents the following
 689         // situation.  Say for example we have a perfect octagon defined by the input winding.
 690         // Say our chopping plane (defined by normal and dist) is essentially the same plane
 691         // that the octagon is sitting on.  Well, due to rounding errors, it may be that point
 692         // 1 of the octagon might be in front, point 2 might be in back, point 3 might be in
 693         // front, point 4 might be in back, and so on.  So we could end up with a very ugly-
 694         // looking chopped winding, and this might be undesirable, and would at least lead to
 695         // a possible exhaustion of MAX_POINTS_ON_WINDING.  It's better to assume that points
 696         // very very close to the plane are on the plane, using an infinitesimal epsilon amount.
 697
 698         // Now, the original ChopWindingInPlace() function used a vec_t-based winding_t.
 699         // So this minimum epsilon is quite similar to casting the higher resolution numbers to
 700         // the lower resolution and comparing them in the lower resolution mode.  We explicitly
 701         // choose the minimum epsilon as something around the vec_t epsilon of one because we
 702         // want the resolution of vec_accu_t to have a large resolution around the epsilon.
 703         // Some of that leftover resolution even goes away after we scale to points far away.
 704
 705         // Here is a further discussion regarding the choice of smallestEpsilonAllowed.
 706         // In the 32 float world (we can assume vec_t is that), the "epsilon around 1.0" is
 707         // 0.00000011921.  In the 64 bit float world (we can assume vec_accu_t is that), the
 708         // "epsilon around 1.0" is 0.00000000000000022204.  (By the way these two epsilons
 709         // are defined as VEC_SMALLEST_EPSILON_AROUND_ONE VEC_ACCU_SMALLEST_EPSILON_AROUND_ONE
 710         // respectively.)  If you divide the first by the second, you get approximately
 711         // 536,885,246.  Dividing that number by 200,000 (a typical base winding coordinate)
 712         // gives 2684.  So in other words, if our smallestEpsilonAllowed was chosen as exactly
 713         // VEC_SMALLEST_EPSILON_AROUND_ONE, you would be guaranteed at least 2000 "ticks" in
 714         // 64-bit land inside of the epsilon for all numbers we're dealing with.
 715
 716         static const vec_accu_t smallestEpsilonAllowed = ( (vec_accu_t) VEC_SMALLEST_EPSILON_AROUND_ONE ) * 0.5;
 717         if ( crudeEpsilon < smallestEpsilonAllowed ) {
 718                 fineEpsilon = smallestEpsilonAllowed;
 719         }
 720         else{fineEpsilon = (vec_accu_t) crudeEpsilon; }
 721
 722         in = *inout;
 723         counts[0] = counts[1] = counts[2] = 0;
 724         VectorCopyRegularToAccu( normal, normalAccu );
 725
 726         for ( i = 0; i < in->numpoints; i++ )
 727         {
 728                 dists[i] = DotProductAccu( in->p[i], normalAccu ) - dist;
 729                 if ( dists[i] > fineEpsilon ) {
 730                         sides[i] = SIDE_FRONT;
 731                 }
 732                 else if ( dists[i] < -fineEpsilon ) {
 733                         sides[i] = SIDE_BACK;
 734                 }
 735                 else{sides[i] = SIDE_ON; }
 736                 counts[sides[i]]++;
 737         }
 738         sides[i] = sides[0];
 739         dists[i] = dists[0];
 740
 741         // I'm wondering if whatever code that handles duplicate planes is robust enough
 742         // that we never get a case where two nearly equal planes result in 2 NULL windings
 743         // due to the 'if' statement below.  TODO: Investigate this.
 744         if ( !counts[SIDE_FRONT] ) {
 745                 FreeWindingAccu( in );
 746                 *inout = NULL;
 747                 return;
 748         }
 749         if ( !counts[SIDE_BACK] ) {
 750                 return; // Winding is unmodified.
 751         }
 752
 753         // NOTE: The least number of points that a winding can have at this point is 2.
 754         // In that case, one point is SIDE_FRONT and the other is SIDE_BACK.
 755
 756         maxpts = counts[SIDE_FRONT] + 2; // We dynamically expand if this is too small.
 757         f = AllocWindingAccu( maxpts );
 758
 759         for ( i = 0; i < in->numpoints; i++ )
 760         {
 761                 p1 = in->p[i];
 762
 763                 if ( sides[i] == SIDE_ON || sides[i] == SIDE_FRONT ) {
 764                         if ( f->numpoints >= MAX_POINTS_ON_WINDING ) {
 765                                 Error( "ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING" );
 766                         }
 767                         if ( f->numpoints >= maxpts ) { // This will probably never happen.
 768                                 Sys_FPrintf( SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n" );
 769                                 f = CopyWindingAccuIncreaseSizeAndFreeOld( f );
 770                                 maxpts++;
 771                         }
 772                         VectorCopyAccu( p1, f->p[f->numpoints] );
 773                         f->numpoints++;
 774                         if ( sides[i] == SIDE_ON ) {
 775                                 continue;
 776                         }
 777                 }
 778                 if ( sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i] ) {
 779                         continue;
 780                 }
 781
 782                 // Generate a split point.
 783                 p2 = in->p[( ( i + 1 ) == in->numpoints ) ? 0 : ( i + 1 )];
 784
 785                 // The divisor's absolute value is greater than the dividend's absolute value.
 786                 // w is in the range (0,1).
 787                 w = dists[i] / ( dists[i] - dists[i + 1] );
 788
 789                 for ( j = 0; j < 3; j++ )
 790                 {
 791                         // Avoid round-off error when possible.  Check axis-aligned normal.
 792                         if ( normal[j] == 1 ) {
 793                                 mid[j] = dist;
 794                         }
 795                         else if ( normal[j] == -1 ) {
 796                                 mid[j] = -dist;
 797                         }
 798                         else{mid[j] = p1[j] + ( w * ( p2[j] - p1[j] ) ); }
 799                 }
 800                 if ( f->numpoints >= MAX_POINTS_ON_WINDING ) {
 801                         Error( "ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING" );
 802                 }
 803                 if ( f->numpoints >= maxpts ) { // This will probably never happen.
 804                         Sys_FPrintf( SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n" );
 805                         f = CopyWindingAccuIncreaseSizeAndFreeOld( f );
 806                         maxpts++;
 807                 }
 808                 VectorCopyAccu( mid, f->p[f->numpoints] );
 809                 f->numpoints++;
 810         }
 811
 812         FreeWindingAccu( in );
 813         *inout = f;
 814 }
 815
 816 /*
 817    =============
 818    ChopWindingInPlace
 819    =============
 820  */
 821 void ChopWindingInPlace( winding_t **inout, vec3_t normal, vec_t dist, vec_t epsilon ){
 822         winding_t   *in;
 823         vec_t dists[MAX_POINTS_ON_WINDING + 4];
 824         int sides[MAX_POINTS_ON_WINDING + 4];
 825         int counts[3];
 826         static vec_t dot;           // VC 4.2 optimizer bug if not static
 827         int i, j;
 828         vec_t   *p1, *p2;
 829         vec3_t mid;
 830         winding_t   *f;
 831         int maxpts;
 832
 833         in = *inout;
 834         counts[0] = counts[1] = counts[2] = 0;
 835
 836 // determine sides for each point
 837         for ( i = 0 ; i < in->numpoints ; i++ )
 838         {
 839                 dot = DotProduct( in->p[i], normal );
 840                 dot -= dist;
 841                 dists[i] = dot;
 842                 if ( dot > epsilon ) {
 843                         sides[i] = SIDE_FRONT;
 844                 }
 845                 else if ( dot < -epsilon ) {
 846                         sides[i] = SIDE_BACK;
 847                 }
 848                 else
 849                 {
 850                         sides[i] = SIDE_ON;
 851                 }
 852                 counts[sides[i]]++;
 853         }
 854         sides[i] = sides[0];
 855         dists[i] = dists[0];
 856
 857         if ( !counts[0] ) {
 858                 FreeWinding( in );
 859                 *inout = NULL;
 860                 return;
 861         }
 862         if ( !counts[1] ) {
 863                 return;     // inout stays the same
 864
 865         }
 866         maxpts = in->numpoints + 4;   // cant use counts[0]+2 because
 867                                       // of fp grouping errors
 868
 869         f = AllocWinding( maxpts );
 870
 871         for ( i = 0 ; i < in->numpoints ; i++ )
 872         {
 873                 p1 = in->p[i];
 874
 875                 if ( sides[i] == SIDE_ON ) {
 876                         VectorCopy( p1, f->p[f->numpoints] );
 877                         f->numpoints++;
 878                         continue;
 879                 }
 880
 881                 if ( sides[i] == SIDE_FRONT ) {
 882                         VectorCopy( p1, f->p[f->numpoints] );
 883                         f->numpoints++;
 884                 }
 885
 886                 if ( sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i] ) {
 887                         continue;
 888                 }
 889
 890                 // generate a split point
 891                 p2 = in->p[( i + 1 ) % in->numpoints];
 892
 893                 dot = dists[i] / ( dists[i] - dists[i + 1] );
 894                 for ( j = 0 ; j < 3 ; j++ )
 895                 {   // avoid round off error when possible
 896                         if ( normal[j] == 1 ) {
 897                                 mid[j] = dist;
 898                         }
 899                         else if ( normal[j] == -1 ) {
 900                                 mid[j] = -dist;
 901                         }
 902                         else{
 903                                 mid[j] = p1[j] + dot * ( p2[j] - p1[j] );
 904                         }
 905                 }
 906
 907                 VectorCopy( mid, f->p[f->numpoints] );
 908                 f->numpoints++;
 909         }
 910
 911         if ( f->numpoints > maxpts ) {
 912                 Error( "ClipWinding: points exceeded estimate" );
 913         }
 914         if ( f->numpoints > MAX_POINTS_ON_WINDING ) {
 915                 Error( "ClipWinding: MAX_POINTS_ON_WINDING" );
 916         }
 917
 918         FreeWinding( in );
 919         *inout = f;
 920 }
 921
 922
 923 /*
 924    =================
 925    ChopWinding
 926
 927    Returns the fragment of in that is on the front side
 928    of the cliping plane.  The original is freed.
 929    =================
 930  */
 931 winding_t   *ChopWinding( winding_t *in, vec3_t normal, vec_t dist ){
 932         winding_t   *f, *b;
 933
 934         ClipWindingEpsilon( in, normal, dist, ON_EPSILON, &f, &b );
 935         FreeWinding( in );
 936         if ( b ) {
 937                 FreeWinding( b );
 938         }
 939         return f;
 940 }
 941
 942
 943 /*
 944    =================
 945    CheckWinding
 946
 947    =================
 948  */
 949 void CheckWinding( winding_t *w ){
 950         int i, j;
 951         vec_t   *p1, *p2;
 952         vec_t d, edgedist;
 953         vec3_t dir, edgenormal, facenormal;
 954         vec_t area;
 955         vec_t facedist;
 956
 957         if ( w->numpoints < 3 ) {
 958                 Error( "CheckWinding: %i points",w->numpoints );
 959         }
 960
 961         area = WindingArea( w );
 962         if ( area < 1 ) {
 963                 Error( "CheckWinding: %f area", area );
 964         }
 965
 966         WindingPlane( w, facenormal, &facedist );
 967
 968         for ( i = 0 ; i < w->numpoints ; i++ )
 969         {
 970                 p1 = w->p[i];
 971
 972                 for ( j = 0 ; j < 3 ; j++ )
 973                         if ( p1[j] > MAX_WORLD_COORD || p1[j] < MIN_WORLD_COORD ) {
 974                                 Error( "CheckFace: MAX_WORLD_COORD exceeded: %f",p1[j] );
 975                         }
 976
 977                 j = i + 1 == w->numpoints ? 0 : i + 1;
 978
 979                 // check the point is on the face plane
 980                 d = DotProduct( p1, facenormal ) - facedist;
 981                 if ( d < -ON_EPSILON || d > ON_EPSILON ) {
 982                         Error( "CheckWinding: point off plane" );
 983                 }
 984
 985                 // check the edge isnt degenerate
 986                 p2 = w->p[j];
 987                 VectorSubtract( p2, p1, dir );
 988
 989                 if ( VectorLength( dir ) < ON_EPSILON ) {
 990                         Error( "CheckWinding: degenerate edge" );
 991                 }
 992
 993                 CrossProduct( facenormal, dir, edgenormal );
 994                 VectorNormalize( edgenormal, edgenormal );
 995                 edgedist = DotProduct( p1, edgenormal );
 996                 edgedist += ON_EPSILON;
 997
 998                 // all other points must be on front side
 999                 for ( j = 0 ; j < w->numpoints ; j++ )
1000                 {
1001                         if ( j == i ) {
1002                                 continue;
1003                         }
1004                         d = DotProduct( w->p[j], edgenormal );
1005                         if ( d > edgedist ) {
1006                                 Error( "CheckWinding: non-convex" );
1007                         }
1008                 }
1009         }
1010 }
1011
1012
1013 /*
1014    ============
1015    WindingOnPlaneSide
1016    ============
1017  */
1018 int     WindingOnPlaneSide( winding_t *w, vec3_t normal, vec_t dist ){
1019         qboolean front, back;
1020         int i;
1021         vec_t d;
1022
1023         front = qfalse;
1024         back = qfalse;
1025         for ( i = 0 ; i < w->numpoints ; i++ )
1026         {
1027                 d = DotProduct( w->p[i], normal ) - dist;
1028                 if ( d < -ON_EPSILON ) {
1029                         if ( front ) {
1030                                 return SIDE_CROSS;
1031                         }
1032                         back = qtrue;
1033                         continue;
1034                 }
1035                 if ( d > ON_EPSILON ) {
1036                         if ( back ) {
1037                                 return SIDE_CROSS;
1038                         }
1039                         front = qtrue;
1040                         continue;
1041                 }
1042         }
1043
1044         if ( back ) {
1045                 return SIDE_BACK;
1046         }
1047         if ( front ) {
1048                 return SIDE_FRONT;
1049         }
1050         return SIDE_ON;
1051 }
1052
1053
1054 /*
1055    =================
1056    AddWindingToConvexHull
1057
1058    Both w and *hull are on the same plane
1059    =================
1060  */
1061 #define MAX_HULL_POINTS     128
1062 void    AddWindingToConvexHull( winding_t *w, winding_t **hull, vec3_t normal ) {
1063         int i, j, k;
1064         float       *p, *copy;
1065         vec3_t dir;
1066         float d;
1067         int numHullPoints, numNew;
1068         vec3_t hullPoints[MAX_HULL_POINTS];
1069         vec3_t newHullPoints[MAX_HULL_POINTS];
1070         vec3_t hullDirs[MAX_HULL_POINTS];
1071         qboolean hullSide[MAX_HULL_POINTS];
1072         qboolean outside;
1073
1074         if ( !*hull ) {
1075                 *hull = CopyWinding( w );
1076                 return;
1077         }
1078
1079         numHullPoints = ( *hull )->numpoints;
1080         memcpy( hullPoints, ( *hull )->p, numHullPoints * sizeof( vec3_t ) );
1081
1082         for ( i = 0 ; i < w->numpoints ; i++ ) {
1083                 p = w->p[i];
1084
1085                 // calculate hull side vectors
1086                 for ( j = 0 ; j < numHullPoints ; j++ ) {
1087                         k = ( j + 1 ) % numHullPoints;
1088
1089                         VectorSubtract( hullPoints[k], hullPoints[j], dir );
1090                         VectorNormalize( dir, dir );
1091                         CrossProduct( normal, dir, hullDirs[j] );
1092                 }
1093
1094                 outside = qfalse;
1095                 for ( j = 0 ; j < numHullPoints ; j++ ) {
1096                         VectorSubtract( p, hullPoints[j], dir );
1097                         d = DotProduct( dir, hullDirs[j] );
1098                         if ( d >= ON_EPSILON ) {
1099                                 outside = qtrue;
1100                         }
1101                         if ( d >= -ON_EPSILON ) {
1102                                 hullSide[j] = qtrue;
1103                         }
1104                         else {
1105                                 hullSide[j] = qfalse;
1106                         }
1107                 }
1108
1109                 // if the point is effectively inside, do nothing
1110                 if ( !outside ) {
1111                         continue;
1112                 }
1113
1114                 // find the back side to front side transition
1115                 for ( j = 0 ; j < numHullPoints ; j++ ) {
1116                         if ( !hullSide[ j % numHullPoints ] && hullSide[ ( j + 1 ) % numHullPoints ] ) {
1117                                 break;
1118                         }
1119                 }
1120                 if ( j == numHullPoints ) {
1121                         continue;
1122                 }
1123
1124                 // insert the point here
1125                 VectorCopy( p, newHullPoints[0] );
1126                 numNew = 1;
1127
1128                 // copy over all points that aren't double fronts
1129                 j = ( j + 1 ) % numHullPoints;
1130                 for ( k = 0 ; k < numHullPoints ; k++ ) {
1131                         if ( hullSide[ ( j + k ) % numHullPoints ] && hullSide[ ( j + k + 1 ) % numHullPoints ] ) {
1132                                 continue;
1133                         }
1134                         copy = hullPoints[ ( j + k + 1 ) % numHullPoints ];
1135                         VectorCopy( copy, newHullPoints[numNew] );
1136                         numNew++;
1137                 }
1138
1139                 numHullPoints = numNew;
1140                 memcpy( hullPoints, newHullPoints, numHullPoints * sizeof( vec3_t ) );
1141         }
1142
1143         FreeWinding( *hull );
1144         w = AllocWinding( numHullPoints );
1145         w->numpoints = numHullPoints;
1146         *hull = w;
1147         memcpy( w->p, hullPoints, numHullPoints * sizeof( vec3_t ) );
1148 }