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