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