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