tools/quake2/extra/common/polylib.c

   1 /*
   2 ===========================================================================
   3 Copyright (C) 1997-2006 Id Software, Inc.
   4
   5 This file is part of Quake 2 Tools source code.
   6
   7 Quake 2 Tools source code is free software; you can redistribute it
   8 and/or modify it under the terms of the GNU General Public License as
   9 published by the Free Software Foundation; either version 2 of the License,
  10 or (at your option) any later version.
  11
  12 Quake 2 Tools source code is distributed in the hope that it will be
  13 useful, 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 Quake 2 Tools source code; 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 "polylib.h"
  26
  27
  28 extern int numthreads;
  29
  30 // counters are only bumped when running single threaded,
  31 // because they are an awefull coherence problem
  32 int     c_active_windings;
  33 int     c_peak_windings;
  34 int     c_winding_allocs;
  35 int     c_winding_points;
  36
  37 #define BOGUS_RANGE     8192
  38
  39 void pw(winding_t *w)
  40 {
  41         int             i;
  42         for (i=0 ; i<w->numpoints ; i++)
  43                 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 {
  54         winding_t       *w;
  55         int                     s;
  56
  57         if (numthreads == 1)
  58         {
  59                 c_winding_allocs++;
  60                 c_winding_points += points;
  61                 c_active_windings++;
  62                 if (c_active_windings > c_peak_windings)
  63                         c_peak_windings = c_active_windings;
  64         }
  65         s = sizeof(vec_t)*3*points + sizeof(int);
  66         w = malloc (s);
  67         memset (w, 0, s);
  68         return w;
  69 }
  70
  71 void FreeWinding (winding_t *w)
  72 {
  73         if (*(unsigned *)w == 0xdeaddead)
  74                 Error ("FreeWinding: freed a freed winding");
  75         *(unsigned *)w = 0xdeaddead;
  76
  77         if (numthreads == 1)
  78                 c_active_windings--;
  79         free (w);
  80 }
  81
  82 /*
  83 ============
  84 RemoveColinearPoints
  85 ============
  86 */
  87 int     c_removed;
  88
  89 void    RemoveColinearPoints (winding_t *w)
  90 {
  91         int             i, j, k;
  92         vec3_t  v1, v2;
  93         int             nump;
  94         vec3_t  p[MAX_POINTS_ON_WINDING];
  95
  96         nump = 0;
  97         for (i=0 ; i<w->numpoints ; i++)
  98         {
  99                 j = (i+1)%w->numpoints;
 100                 k = (i+w->numpoints-1)%w->numpoints;
 101                 VectorSubtract (w->p[j], w->p[i], v1);
 102                 VectorSubtract (w->p[i], w->p[k], v2);
 103                 VectorNormalize(v1,v1);
 104                 VectorNormalize(v2,v2);
 105                 if (DotProduct(v1, v2) < 0.999)
 106                 {
 107                         VectorCopy (w->p[i], p[nump]);
 108                         nump++;
 109                 }
 110         }
 111
 112         if (nump == w->numpoints)
 113                 return;
 114
 115         if (numthreads == 1)
 116                 c_removed += w->numpoints - nump;
 117         w->numpoints = nump;
 118         memcpy (w->p, p, nump*sizeof(p[0]));
 119 }
 120
 121 /*
 122 ============
 123 WindingPlane
 124 ============
 125 */
 126 void WindingPlane (winding_t *w, vec3_t normal, vec_t *dist)
 127 {
 128         vec3_t  v1, v2;
 129
 130         VectorSubtract (w->p[1], w->p[0], v1);
 131         VectorSubtract (w->p[2], w->p[0], v2);
 132         CrossProduct (v2, v1, normal);
 133         VectorNormalize (normal, normal);
 134         *dist = DotProduct (w->p[0], normal);
 135
 136 }
 137
 138 /*
 139 =============
 140 WindingArea
 141 =============
 142 */
 143 vec_t   WindingArea (winding_t *w)
 144 {
 145         int             i;
 146         vec3_t  d1, d2, cross;
 147         vec_t   total;
 148
 149         total = 0;
 150         for (i=2 ; i<w->numpoints ; i++)
 151         {
 152                 VectorSubtract (w->p[i-1], w->p[0], d1);
 153                 VectorSubtract (w->p[i], w->p[0], d2);
 154                 CrossProduct (d1, d2, cross);
 155                 total += 0.5 * VectorLength ( cross );
 156         }
 157         return total;
 158 }
 159
 160 void    WindingBounds (winding_t *w, vec3_t mins, vec3_t maxs)
 161 {
 162         vec_t   v;
 163         int             i,j;
 164
 165         mins[0] = mins[1] = mins[2] = 99999;
 166         maxs[0] = maxs[1] = maxs[2] = -99999;
 167
 168         for (i=0 ; i<w->numpoints ; i++)
 169         {
 170                 for (j=0 ; j<3 ; j++)
 171                 {
 172                         v = w->p[i][j];
 173                         if (v < mins[j])
 174                                 mins[j] = v;
 175                         if (v > maxs[j])
 176                                 maxs[j] = v;
 177                 }
 178         }
 179 }
 180
 181 /*
 182 =============
 183 WindingCenter
 184 =============
 185 */
 186 void    WindingCenter (winding_t *w, vec3_t center)
 187 {
 188         int             i;
 189         float   scale;
 190
 191         VectorCopy (vec3_origin, center);
 192         for (i=0 ; i<w->numpoints ; i++)
 193                 VectorAdd (w->p[i], center, center);
 194
 195         scale = 1.0/w->numpoints;
 196         VectorScale (center, scale, center);
 197 }
 198
 199 /*
 200 =================
 201 BaseWindingForPlane
 202 =================
 203 */
 204 winding_t *BaseWindingForPlane (vec3_t normal, vec_t dist)
 205 {
 206         int             i, x;
 207         vec_t   max, v;
 208         vec3_t  org, vright, vup;
 209         winding_t       *w;
 210
 211 // find the major axis
 212
 213         max = -BOGUS_RANGE;
 214         x = -1;
 215         for (i=0 ; i<3; i++)
 216         {
 217                 v = fabs(normal[i]);
 218                 if (v > max)
 219                 {
 220                         x = i;
 221                         max = v;
 222                 }
 223         }
 224         if (x==-1)
 225                 Error ("BaseWindingForPlane: no axis found");
 226
 227         VectorCopy (vec3_origin, vup);
 228         switch (x)
 229         {
 230         case 0:
 231         case 1:
 232                 vup[2] = 1;
 233                 break;
 234         case 2:
 235                 vup[0] = 1;
 236                 break;
 237         }
 238
 239         v = DotProduct (vup, normal);
 240         VectorMA (vup, -v, normal, vup);
 241         VectorNormalize (vup, vup);
 242
 243         VectorScale (normal, dist, org);
 244
 245         CrossProduct (vup, normal, vright);
 246
 247         VectorScale (vup, 8192, vup);
 248         VectorScale (vright, 8192, vright);
 249
 250 // project a really big axis aligned box onto the plane
 251         w = AllocWinding (4);
 252
 253         VectorSubtract (org, vright, w->p[0]);
 254         VectorAdd (w->p[0], vup, w->p[0]);
 255
 256         VectorAdd (org, vright, w->p[1]);
 257         VectorAdd (w->p[1], vup, w->p[1]);
 258
 259         VectorAdd (org, vright, w->p[2]);
 260         VectorSubtract (w->p[2], vup, w->p[2]);
 261
 262         VectorSubtract (org, vright, w->p[3]);
 263         VectorSubtract (w->p[3], vup, w->p[3]);
 264
 265         w->numpoints = 4;
 266
 267         return w;
 268 }
 269
 270 /*
 271 ==================
 272 CopyWinding
 273 ==================
 274 */
 275 winding_t       *CopyWinding (winding_t *w)
 276 {
 277         int                     size;
 278         winding_t       *c;
 279
 280         c = AllocWinding (w->numpoints);
 281         size = (int)((winding_t *)0)->p[w->numpoints];
 282         memcpy (c, w, size);
 283         return c;
 284 }
 285
 286 /*
 287 ==================
 288 ReverseWinding
 289 ==================
 290 */
 291 winding_t       *ReverseWinding (winding_t *w)
 292 {
 293         int                     i;
 294         winding_t       *c;
 295
 296         c = AllocWinding (w->numpoints);
 297         for (i=0 ; i<w->numpoints ; i++)
 298         {
 299                 VectorCopy (w->p[w->numpoints-1-i], c->p[i]);
 300         }
 301         c->numpoints = w->numpoints;
 302         return c;
 303 }
 304
 305
 306 /*
 307 =============
 308 ClipWindingEpsilon
 309 =============
 310 */
 311 void    ClipWindingEpsilon (winding_t *in, vec3_t normal, vec_t dist,
 312                                 vec_t epsilon, winding_t **front, winding_t **back)
 313 {
 314         vec_t   dists[MAX_POINTS_ON_WINDING+4];
 315         int             sides[MAX_POINTS_ON_WINDING+4];
 316         int             counts[3];
 317         static  vec_t   dot;            // VC 4.2 optimizer bug if not static
 318         int             i, j;
 319         vec_t   *p1, *p2;
 320         vec3_t  mid;
 321         winding_t       *f, *b;
 322         int             maxpts;
 323
 324         counts[0] = counts[1] = counts[2] = 0;
 325
 326 // determine sides for each point
 327         for (i=0 ; i<in->numpoints ; i++)
 328         {
 329                 dot = DotProduct (in->p[i], normal);
 330                 dot -= dist;
 331                 dists[i] = dot;
 332                 if (dot > epsilon)
 333                         sides[i] = SIDE_FRONT;
 334                 else if (dot < -epsilon)
 335                         sides[i] = SIDE_BACK;
 336                 else
 337                 {
 338                         sides[i] = SIDE_ON;
 339                 }
 340                 counts[sides[i]]++;
 341         }
 342         sides[i] = sides[0];
 343         dists[i] = dists[0];
 344
 345         *front = *back = NULL;
 346
 347         if (!counts[0])
 348         {
 349                 *back = CopyWinding (in);
 350                 return;
 351         }
 352         if (!counts[1])
 353         {
 354                 *front = CopyWinding (in);
 355                 return;
 356         }
 357
 358         maxpts = in->numpoints+4;       // cant use counts[0]+2 because
 359                                                                 // of fp grouping errors
 360
 361         *front = f = AllocWinding (maxpts);
 362         *back = b = AllocWinding (maxpts);
 363
 364         for (i=0 ; i<in->numpoints ; i++)
 365         {
 366                 p1 = in->p[i];
 367
 368                 if (sides[i] == SIDE_ON)
 369                 {
 370                         VectorCopy (p1, f->p[f->numpoints]);
 371                         f->numpoints++;
 372                         VectorCopy (p1, b->p[b->numpoints]);
 373                         b->numpoints++;
 374                         continue;
 375                 }
 376
 377                 if (sides[i] == SIDE_FRONT)
 378                 {
 379                         VectorCopy (p1, f->p[f->numpoints]);
 380                         f->numpoints++;
 381                 }
 382                 if (sides[i] == SIDE_BACK)
 383                 {
 384                         VectorCopy (p1, b->p[b->numpoints]);
 385                         b->numpoints++;
 386                 }
 387
 388                 if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
 389                         continue;
 390
 391         // generate a split point
 392                 p2 = in->p[(i+1)%in->numpoints];
 393
 394                 dot = dists[i] / (dists[i]-dists[i+1]);
 395                 for (j=0 ; j<3 ; j++)
 396                 {       // avoid round off error when possible
 397                         if (normal[j] == 1)
 398                                 mid[j] = dist;
 399                         else if (normal[j] == -1)
 400                                 mid[j] = -dist;
 401                         else
 402                                 mid[j] = p1[j] + dot*(p2[j]-p1[j]);
 403                 }
 404
 405                 VectorCopy (mid, f->p[f->numpoints]);
 406                 f->numpoints++;
 407                 VectorCopy (mid, b->p[b->numpoints]);
 408                 b->numpoints++;
 409         }
 410
 411         if (f->numpoints > maxpts || b->numpoints > maxpts)
 412                 Error ("ClipWinding: points exceeded estimate");
 413         if (f->numpoints > MAX_POINTS_ON_WINDING || b->numpoints > MAX_POINTS_ON_WINDING)
 414                 Error ("ClipWinding: MAX_POINTS_ON_WINDING");
 415 }
 416
 417
 418 /*
 419 =============
 420 ChopWindingInPlace
 421 =============
 422 */
 423 void ChopWindingInPlace (winding_t **inout, vec3_t normal, vec_t dist, vec_t epsilon)
 424 {
 425         winding_t       *in;
 426         vec_t   dists[MAX_POINTS_ON_WINDING+4];
 427         int             sides[MAX_POINTS_ON_WINDING+4];
 428         int             counts[3];
 429         static  vec_t   dot;            // VC 4.2 optimizer bug if not static
 430         int             i, j;
 431         vec_t   *p1, *p2;
 432         vec3_t  mid;
 433         winding_t       *f;
 434         int             maxpts;
 435
 436         in = *inout;
 437         counts[0] = counts[1] = counts[2] = 0;
 438
 439 // determine sides for each point
 440         for (i=0 ; i<in->numpoints ; i++)
 441         {
 442                 dot = DotProduct (in->p[i], normal);
 443                 dot -= dist;
 444                 dists[i] = dot;
 445                 if (dot > epsilon)
 446                         sides[i] = SIDE_FRONT;
 447                 else if (dot < -epsilon)
 448                         sides[i] = SIDE_BACK;
 449                 else
 450                 {
 451                         sides[i] = SIDE_ON;
 452                 }
 453                 counts[sides[i]]++;
 454         }
 455         sides[i] = sides[0];
 456         dists[i] = dists[0];
 457
 458         if (!counts[0])
 459         {
 460                 FreeWinding (in);
 461                 *inout = NULL;
 462                 return;
 463         }
 464         if (!counts[1])
 465                 return;         // inout stays the same
 466
 467         maxpts = in->numpoints+4;       // cant use counts[0]+2 because
 468                                                                 // of fp grouping errors
 469
 470         f = AllocWinding (maxpts);
 471
 472         for (i=0 ; i<in->numpoints ; i++)
 473         {
 474                 p1 = in->p[i];
 475
 476                 if (sides[i] == SIDE_ON)
 477                 {
 478                         VectorCopy (p1, f->p[f->numpoints]);
 479                         f->numpoints++;
 480                         continue;
 481                 }
 482
 483                 if (sides[i] == SIDE_FRONT)
 484                 {
 485                         VectorCopy (p1, f->p[f->numpoints]);
 486                         f->numpoints++;
 487                 }
 488
 489                 if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
 490                         continue;
 491
 492         // generate a split point
 493                 p2 = in->p[(i+1)%in->numpoints];
 494
 495                 dot = dists[i] / (dists[i]-dists[i+1]);
 496                 for (j=0 ; j<3 ; j++)
 497                 {       // avoid round off error when possible
 498                         if (normal[j] == 1)
 499                                 mid[j] = dist;
 500                         else if (normal[j] == -1)
 501                                 mid[j] = -dist;
 502                         else
 503                                 mid[j] = p1[j] + dot*(p2[j]-p1[j]);
 504                 }
 505
 506                 VectorCopy (mid, f->p[f->numpoints]);
 507                 f->numpoints++;
 508         }
 509
 510         if (f->numpoints > maxpts)
 511                 Error ("ClipWinding: points exceeded estimate");
 512         if (f->numpoints > MAX_POINTS_ON_WINDING)
 513                 Error ("ClipWinding: MAX_POINTS_ON_WINDING");
 514
 515         FreeWinding (in);
 516         *inout = f;
 517 }
 518
 519
 520 /*
 521 =================
 522 ChopWinding
 523
 524 Returns the fragment of in that is on the front side
 525 of the cliping plane.  The original is freed.
 526 =================
 527 */
 528 winding_t       *ChopWinding (winding_t *in, vec3_t normal, vec_t dist)
 529 {
 530         winding_t       *f, *b;
 531
 532         ClipWindingEpsilon (in, normal, dist, ON_EPSILON, &f, &b);
 533         FreeWinding (in);
 534         if (b)
 535                 FreeWinding (b);
 536         return f;
 537 }
 538
 539
 540 /*
 541 =================
 542 CheckWinding
 543
 544 =================
 545 */
 546 void CheckWinding (winding_t *w)
 547 {
 548         int             i, j;
 549         vec_t   *p1, *p2;
 550         vec_t   d, edgedist;
 551         vec3_t  dir, edgenormal, facenormal;
 552         vec_t   area;
 553         vec_t   facedist;
 554
 555         if (w->numpoints < 3)
 556                 Error ("CheckWinding: %i points",w->numpoints);
 557
 558         area = WindingArea(w);
 559         if (area < 1)
 560                 Error ("CheckWinding: %f area", area);
 561
 562         WindingPlane (w, facenormal, &facedist);
 563
 564         for (i=0 ; i<w->numpoints ; i++)
 565         {
 566                 p1 = w->p[i];
 567
 568                 for (j=0 ; j<3 ; j++)
 569                         if (p1[j] > BOGUS_RANGE || p1[j] < -BOGUS_RANGE)
 570                                 Error ("CheckFace: BUGUS_RANGE: %f",p1[j]);
 571
 572                 j = i+1 == w->numpoints ? 0 : i+1;
 573
 574         // check the point is on the face plane
 575                 d = DotProduct (p1, facenormal) - facedist;
 576                 if (d < -ON_EPSILON || d > ON_EPSILON)
 577                         Error ("CheckWinding: point off plane");
 578
 579         // check the edge isnt degenerate
 580                 p2 = w->p[j];
 581                 VectorSubtract (p2, p1, dir);
 582
 583                 if (VectorLength (dir) < ON_EPSILON)
 584                         Error ("CheckWinding: degenerate edge");
 585
 586                 CrossProduct (facenormal, dir, edgenormal);
 587                 VectorNormalize (edgenormal, edgenormal);
 588                 edgedist = DotProduct (p1, edgenormal);
 589                 edgedist += ON_EPSILON;
 590
 591         // all other points must be on front side
 592                 for (j=0 ; j<w->numpoints ; j++)
 593                 {
 594                         if (j == i)
 595                                 continue;
 596                         d = DotProduct (w->p[j], edgenormal);
 597                         if (d > edgedist)
 598                                 Error ("CheckWinding: non-convex");
 599                 }
 600         }
 601 }
 602
 603
 604 /*
 605 ============
 606 WindingOnPlaneSide
 607 ============
 608 */
 609 int             WindingOnPlaneSide (winding_t *w, vec3_t normal, vec_t dist)
 610 {
 611         qboolean        front, back;
 612         int                     i;
 613         vec_t           d;
 614
 615         front = false;
 616         back = false;
 617         for (i=0 ; i<w->numpoints ; i++)
 618         {
 619                 d = DotProduct (w->p[i], normal) - dist;
 620                 if (d < -ON_EPSILON)
 621                 {
 622                         if (front)
 623                                 return SIDE_CROSS;
 624                         back = true;
 625                         continue;
 626                 }
 627                 if (d > ON_EPSILON)
 628                 {
 629                         if (back)
 630                                 return SIDE_CROSS;
 631                         front = true;
 632                         continue;
 633                 }
 634         }
 635
 636         if (back)
 637                 return SIDE_BACK;
 638         if (front)
 639                 return SIDE_FRONT;
 640         return SIDE_ON;
 641 }
 642