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 void FreeWinding (winding_t *w)
  77 {
  78         if (*(unsigned *)w == 0xdeaddead)
  79                 Error ("FreeWinding: freed a freed winding");
  80         *(unsigned *)w = 0xdeaddead;
  81
  82         if (numthreads == 1)
  83                 c_active_windings--;
  84         free (w);
  85 }
  86
  87 /*
  88 ============
  89 RemoveColinearPoints
  90 ============
  91 */
  92 int     c_removed;
  93
  94 void    RemoveColinearPoints (winding_t *w)
  95 {
  96         int             i, j, k;
  97         vec3_t  v1, v2;
  98         int             nump;
  99         vec3_t  p[MAX_POINTS_ON_WINDING];
 100
 101         nump = 0;
 102         for (i=0 ; i<w->numpoints ; i++)
 103         {
 104                 j = (i+1)%w->numpoints;
 105                 k = (i+w->numpoints-1)%w->numpoints;
 106                 VectorSubtract (w->p[j], w->p[i], v1);
 107                 VectorSubtract (w->p[i], w->p[k], v2);
 108                 VectorNormalize(v1,v1);
 109                 VectorNormalize(v2,v2);
 110                 if (DotProduct(v1, v2) < 0.999)
 111                 {
 112                         VectorCopy (w->p[i], p[nump]);
 113                         nump++;
 114                 }
 115         }
 116
 117         if (nump == w->numpoints)
 118                 return;
 119
 120         if (numthreads == 1)
 121                 c_removed += w->numpoints - nump;
 122         w->numpoints = nump;
 123         memcpy (w->p, p, nump*sizeof(p[0]));
 124 }
 125
 126 /*
 127 ============
 128 WindingPlane
 129 ============
 130 */
 131 void WindingPlane (winding_t *w, vec3_t normal, vec_t *dist)
 132 {
 133         vec3_t  v1, v2;
 134
 135         VectorSubtract (w->p[1], w->p[0], v1);
 136         VectorSubtract (w->p[2], w->p[0], v2);
 137         CrossProduct (v2, v1, normal);
 138         VectorNormalize (normal, normal);
 139         *dist = DotProduct (w->p[0], normal);
 140
 141 }
 142
 143 /*
 144 =============
 145 WindingArea
 146 =============
 147 */
 148 vec_t   WindingArea (winding_t *w)
 149 {
 150         int             i;
 151         vec3_t  d1, d2, cross;
 152         vec_t   total;
 153
 154         total = 0;
 155         for (i=2 ; i<w->numpoints ; i++)
 156         {
 157                 VectorSubtract (w->p[i-1], w->p[0], d1);
 158                 VectorSubtract (w->p[i], w->p[0], d2);
 159                 CrossProduct (d1, d2, cross);
 160                 total += 0.5 * VectorLength ( cross );
 161         }
 162         return total;
 163 }
 164
 165 void    WindingBounds (winding_t *w, vec3_t mins, vec3_t maxs)
 166 {
 167         vec_t   v;
 168         int             i,j;
 169
 170         mins[0] = mins[1] = mins[2] = 99999;
 171         maxs[0] = maxs[1] = maxs[2] = -99999;
 172
 173         for (i=0 ; i<w->numpoints ; i++)
 174         {
 175                 for (j=0 ; j<3 ; j++)
 176                 {
 177                         v = w->p[i][j];
 178                         if (v < mins[j])
 179                                 mins[j] = v;
 180                         if (v > maxs[j])
 181                                 maxs[j] = v;
 182                 }
 183         }
 184 }
 185
 186 /*
 187 =============
 188 WindingCenter
 189 =============
 190 */
 191 void    WindingCenter (winding_t *w, vec3_t center)
 192 {
 193         int             i;
 194         float   scale;
 195
 196         VectorCopy (vec3_origin, center);
 197         for (i=0 ; i<w->numpoints ; i++)
 198                 VectorAdd (w->p[i], center, center);
 199
 200         scale = 1.0/w->numpoints;
 201         VectorScale (center, scale, center);
 202 }
 203
 204 /*
 205 =================
 206 BaseWindingForPlane
 207 =================
 208 */
 209 winding_t *BaseWindingForPlane (vec3_t normal, vec_t dist)
 210 {
 211         int             i, x;
 212         vec_t   max, v;
 213         vec3_t  org, vright, vup;
 214         winding_t       *w;
 215
 216 // find the major axis
 217
 218         max = -BOGUS_RANGE;
 219         x = -1;
 220         for (i=0 ; i<3; i++)
 221         {
 222                 v = fabs(normal[i]);
 223                 if (v > max)
 224                 {
 225                         x = i;
 226                         max = v;
 227                 }
 228         }
 229         if (x==-1)
 230                 Error ("BaseWindingForPlane: no axis found");
 231
 232         VectorCopy (vec3_origin, vup);
 233         switch (x)
 234         {
 235         case 0:
 236         case 1:
 237                 vup[2] = 1;
 238                 break;
 239         case 2:
 240                 vup[0] = 1;
 241                 break;
 242         }
 243
 244         v = DotProduct (vup, normal);
 245         VectorMA (vup, -v, normal, vup);
 246         VectorNormalize (vup, vup);
 247
 248         VectorScale (normal, dist, org);
 249
 250         CrossProduct (vup, normal, vright);
 251
 252         // LordHavoc: this has to use *2 because otherwise some created points may
 253         // be inside the world (think of a diagonal case), and any brush with such
 254         // points should be removed, failure to detect such cases is disasterous
 255         VectorScale (vup, MAX_WORLD_COORD*2, vup);
 256         VectorScale (vright, MAX_WORLD_COORD*2, vright);
 257
 258   // project a really big       axis aligned box onto the plane
 259         w = AllocWinding (4);
 260
 261         VectorSubtract (org, vright, w->p[0]);
 262         VectorAdd (w->p[0], vup, w->p[0]);
 263
 264         VectorAdd (org, vright, w->p[1]);
 265         VectorAdd (w->p[1], vup, w->p[1]);
 266
 267         VectorAdd (org, vright, w->p[2]);
 268         VectorSubtract (w->p[2], vup, w->p[2]);
 269
 270         VectorSubtract (org, vright, w->p[3]);
 271         VectorSubtract (w->p[3], vup, w->p[3]);
 272
 273         w->numpoints = 4;
 274
 275         return w;
 276 }
 277
 278 /*
 279 ==================
 280 CopyWinding
 281 ==================
 282 */
 283 winding_t       *CopyWinding (winding_t *w)
 284 {
 285         int                     size;
 286         winding_t       *c;
 287
 288         c = AllocWinding (w->numpoints);
 289         size = (int)((winding_t *)0)->p[w->numpoints];
 290         memcpy (c, w, size);
 291         return c;
 292 }
 293
 294 /*
 295 ==================
 296 ReverseWinding
 297 ==================
 298 */
 299 winding_t       *ReverseWinding (winding_t *w)
 300 {
 301         int                     i;
 302         winding_t       *c;
 303
 304         c = AllocWinding (w->numpoints);
 305         for (i=0 ; i<w->numpoints ; i++)
 306         {
 307                 VectorCopy (w->p[w->numpoints-1-i], c->p[i]);
 308         }
 309         c->numpoints = w->numpoints;
 310         return c;
 311 }
 312
 313
 314 /*
 315 =============
 316 ClipWindingEpsilon
 317 =============
 318 */
 319 void    ClipWindingEpsilon (winding_t *in, vec3_t normal, vec_t dist,
 320                                 vec_t epsilon, winding_t **front, winding_t **back)
 321 {
 322         vec_t   dists[MAX_POINTS_ON_WINDING+4];
 323         int             sides[MAX_POINTS_ON_WINDING+4];
 324         int             counts[3];
 325         static  vec_t   dot;            // VC 4.2 optimizer bug if not static
 326         int             i, j;
 327         vec_t   *p1, *p2;
 328         vec3_t  mid;
 329         winding_t       *f, *b;
 330         int             maxpts;
 331
 332         counts[0] = counts[1] = counts[2] = 0;
 333
 334 // determine sides for each point
 335         for (i=0 ; i<in->numpoints ; i++)
 336   {
 337
 338                 dot = DotProduct (in->p[i], normal);
 339                 dot -= dist;
 340                 dists[i] = dot;
 341                 if (dot > epsilon)
 342                         sides[i] = SIDE_FRONT;
 343                 else if (dot < -epsilon)
 344                         sides[i] = SIDE_BACK;
 345                 else
 346                 {
 347                         sides[i] = SIDE_ON;
 348                 }
 349                 counts[sides[i]]++;
 350         }
 351         sides[i] = sides[0];
 352         dists[i] = dists[0];
 353
 354         *front = *back = NULL;
 355
 356         if (!counts[0])
 357         {
 358                 *back = CopyWinding (in);
 359                 return;
 360         }
 361         if (!counts[1])
 362         {
 363                 *front = CopyWinding (in);
 364                 return;
 365         }
 366
 367         maxpts = in->numpoints+4;       // cant use counts[0]+2 because
 368                                                                 // of fp grouping errors
 369
 370         *front = f = AllocWinding (maxpts);
 371         *back = b = AllocWinding (maxpts);
 372
 373         for (i=0 ; i<in->numpoints ; i++)
 374         {
 375                 p1 = in->p[i];
 376
 377                 if (sides[i] == SIDE_ON)
 378                 {
 379                         VectorCopy (p1, f->p[f->numpoints]);
 380                         f->numpoints++;
 381                         VectorCopy (p1, b->p[b->numpoints]);
 382                         b->numpoints++;
 383                         continue;
 384                 }
 385
 386                 if (sides[i] == SIDE_FRONT)
 387                 {
 388                         VectorCopy (p1, f->p[f->numpoints]);
 389                         f->numpoints++;
 390                 }
 391                 if (sides[i] == SIDE_BACK)
 392                 {
 393                         VectorCopy (p1, b->p[b->numpoints]);
 394                         b->numpoints++;
 395                 }
 396
 397                 if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
 398                         continue;
 399
 400         // generate a split point
 401                 p2 = in->p[(i+1)%in->numpoints];
 402
 403                 dot = dists[i] / (dists[i]-dists[i+1]);
 404                 for (j=0 ; j<3 ; j++)
 405                 {       // avoid round off error when possible
 406                         if (normal[j] == 1)
 407                                 mid[j] = dist;
 408                         else if (normal[j] == -1)
 409                                 mid[j] = -dist;
 410                         else
 411                                 mid[j] = p1[j] + dot*(p2[j]-p1[j]);
 412                 }
 413
 414                 VectorCopy (mid, f->p[f->numpoints]);
 415                 f->numpoints++;
 416                 VectorCopy (mid, b->p[b->numpoints]);
 417                 b->numpoints++;
 418         }
 419
 420         if (f->numpoints > maxpts || b->numpoints > maxpts)
 421                 Error ("ClipWinding: points exceeded estimate");
 422         if (f->numpoints > MAX_POINTS_ON_WINDING || b->numpoints > MAX_POINTS_ON_WINDING)
 423                 Error ("ClipWinding: MAX_POINTS_ON_WINDING");
 424 }
 425
 426
 427 /*
 428 =============
 429 ChopWindingInPlace
 430 =============
 431 */
 432 void ChopWindingInPlace (winding_t **inout, vec3_t normal, vec_t dist, vec_t epsilon)
 433 {
 434         winding_t       *in;
 435         vec_t   dists[MAX_POINTS_ON_WINDING+4];
 436         int             sides[MAX_POINTS_ON_WINDING+4];
 437         int             counts[3];
 438         static  vec_t   dot;            // VC 4.2 optimizer bug if not static
 439         int             i, j;
 440         vec_t   *p1, *p2;
 441         vec3_t  mid;
 442         winding_t       *f;
 443         int             maxpts;
 444
 445         in = *inout;
 446         counts[0] = counts[1] = counts[2] = 0;
 447
 448 // determine sides for each point
 449         for (i=0 ; i<in->numpoints ; i++)
 450         {
 451                 dot = DotProduct (in->p[i], normal);
 452                 dot -= dist;
 453                 dists[i] = dot;
 454                 if (dot > epsilon)
 455                         sides[i] = SIDE_FRONT;
 456                 else if (dot < -epsilon)
 457                         sides[i] = SIDE_BACK;
 458                 else
 459                 {
 460                         sides[i] = SIDE_ON;
 461                 }
 462                 counts[sides[i]]++;
 463         }
 464         sides[i] = sides[0];
 465         dists[i] = dists[0];
 466
 467         if (!counts[0])
 468         {
 469                 FreeWinding (in);
 470                 *inout = NULL;
 471                 return;
 472         }
 473         if (!counts[1])
 474                 return;         // inout stays the same
 475
 476         maxpts = in->numpoints+4;       // cant use counts[0]+2 because
 477                                                                 // of fp grouping errors
 478
 479         f = AllocWinding (maxpts);
 480
 481         for (i=0 ; i<in->numpoints ; i++)
 482         {
 483                 p1 = in->p[i];
 484
 485                 if (sides[i] == SIDE_ON)
 486                 {
 487                         VectorCopy (p1, f->p[f->numpoints]);
 488                         f->numpoints++;
 489                         continue;
 490                 }
 491
 492                 if (sides[i] == SIDE_FRONT)
 493                 {
 494                         VectorCopy (p1, f->p[f->numpoints]);
 495                         f->numpoints++;
 496                 }
 497
 498                 if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
 499                         continue;
 500
 501         // generate a split point
 502                 p2 = in->p[(i+1)%in->numpoints];
 503
 504                 dot = dists[i] / (dists[i]-dists[i+1]);
 505                 for (j=0 ; j<3 ; j++)
 506                 {       // avoid round off error when possible
 507                         if (normal[j] == 1)
 508                                 mid[j] = dist;
 509                         else if (normal[j] == -1)
 510                                 mid[j] = -dist;
 511                         else
 512                                 mid[j] = p1[j] + dot*(p2[j]-p1[j]);
 513                 }
 514
 515                 VectorCopy (mid, f->p[f->numpoints]);
 516                 f->numpoints++;
 517         }
 518
 519         if (f->numpoints > maxpts)
 520                 Error ("ClipWinding: points exceeded estimate");
 521         if (f->numpoints > MAX_POINTS_ON_WINDING)
 522                 Error ("ClipWinding: MAX_POINTS_ON_WINDING");
 523
 524         FreeWinding (in);
 525         *inout = f;
 526 }
 527
 528
 529 /*
 530 =================
 531 ChopWinding
 532
 533 Returns the fragment of in that is on the front side
 534 of the cliping plane.  The original is freed.
 535 =================
 536 */
 537 winding_t       *ChopWinding (winding_t *in, vec3_t normal, vec_t dist)
 538 {
 539         winding_t       *f, *b;
 540
 541         ClipWindingEpsilon (in, normal, dist, ON_EPSILON, &f, &b);
 542         FreeWinding (in);
 543         if (b)
 544                 FreeWinding (b);
 545         return f;
 546 }
 547
 548
 549 /*
 550 =================
 551 CheckWinding
 552
 553 =================
 554 */
 555 void CheckWinding (winding_t *w)
 556 {
 557         int             i, j;
 558         vec_t   *p1, *p2;
 559         vec_t   d, edgedist;
 560         vec3_t  dir, edgenormal, facenormal;
 561         vec_t   area;
 562         vec_t   facedist;
 563
 564         if (w->numpoints < 3)
 565                 Error ("CheckWinding: %i points",w->numpoints);
 566
 567         area = WindingArea(w);
 568         if (area < 1)
 569                 Error ("CheckWinding: %f area", area);
 570
 571         WindingPlane (w, facenormal, &facedist);
 572
 573         for (i=0 ; i<w->numpoints ; i++)
 574         {
 575                 p1 = w->p[i];
 576
 577                 for (j=0 ; j<3 ; j++)
 578                         if (p1[j] > MAX_WORLD_COORD || p1[j] < MIN_WORLD_COORD)
 579                                 Error ("CheckFace: MAX_WORLD_COORD exceeded: %f",p1[j]);
 580
 581                 j = i+1 == w->numpoints ? 0 : i+1;
 582
 583         // check the point is on the face plane
 584                 d = DotProduct (p1, facenormal) - facedist;
 585                 if (d < -ON_EPSILON || d > ON_EPSILON)
 586                         Error ("CheckWinding: point off plane");
 587
 588         // check the edge isnt degenerate
 589                 p2 = w->p[j];
 590                 VectorSubtract (p2, p1, dir);
 591
 592                 if (VectorLength (dir) < ON_EPSILON)
 593                         Error ("CheckWinding: degenerate edge");
 594
 595                 CrossProduct (facenormal, dir, edgenormal);
 596                 VectorNormalize (edgenormal, edgenormal);
 597                 edgedist = DotProduct (p1, edgenormal);
 598                 edgedist += ON_EPSILON;
 599
 600         // all other points must be on front side
 601                 for (j=0 ; j<w->numpoints ; j++)
 602                 {
 603                         if (j == i)
 604                                 continue;
 605                         d = DotProduct (w->p[j], edgenormal);
 606                         if (d > edgedist)
 607                                 Error ("CheckWinding: non-convex");
 608                 }
 609         }
 610 }
 611
 612
 613 /*
 614 ============
 615 WindingOnPlaneSide
 616 ============
 617 */
 618 int             WindingOnPlaneSide (winding_t *w, vec3_t normal, vec_t dist)
 619 {
 620         qboolean        front, back;
 621         int                     i;
 622         vec_t           d;
 623
 624         front = qfalse;
 625         back = qfalse;
 626         for (i=0 ; i<w->numpoints ; i++)
 627         {
 628                 d = DotProduct (w->p[i], normal) - dist;
 629                 if (d < -ON_EPSILON)
 630                 {
 631                         if (front)
 632                                 return SIDE_CROSS;
 633                         back = qtrue;
 634                         continue;
 635                 }
 636                 if (d > ON_EPSILON)
 637                 {
 638                         if (back)
 639                                 return SIDE_CROSS;
 640                         front = qtrue;
 641                         continue;
 642                 }
 643         }
 644
 645         if (back)
 646                 return SIDE_BACK;
 647         if (front)
 648                 return SIDE_FRONT;
 649         return SIDE_ON;
 650 }
 651
 652
 653 /*
 654 =================
 655 AddWindingToConvexHull
 656
 657 Both w and *hull are on the same plane
 658 =================
 659 */
 660 #define MAX_HULL_POINTS         128
 661 void    AddWindingToConvexHull( winding_t *w, winding_t **hull, vec3_t normal ) {
 662         int                     i, j, k;
 663         float           *p, *copy;
 664         vec3_t          dir;
 665         float           d;
 666         int                     numHullPoints, numNew;
 667         vec3_t          hullPoints[MAX_HULL_POINTS];
 668         vec3_t          newHullPoints[MAX_HULL_POINTS];
 669         vec3_t          hullDirs[MAX_HULL_POINTS];
 670         qboolean        hullSide[MAX_HULL_POINTS];
 671         qboolean        outside;
 672
 673         if ( !*hull ) {
 674                 *hull = CopyWinding( w );
 675                 return;
 676         }
 677
 678         numHullPoints = (*hull)->numpoints;
 679         memcpy( hullPoints, (*hull)->p, numHullPoints * sizeof(vec3_t) );
 680
 681         for ( i = 0 ; i < w->numpoints ; i++ ) {
 682                 p = w->p[i];
 683
 684                 // calculate hull side vectors
 685                 for ( j = 0 ; j < numHullPoints ; j++ ) {
 686                         k = ( j + 1 ) % numHullPoints;
 687
 688                         VectorSubtract( hullPoints[k], hullPoints[j], dir );
 689                         VectorNormalize( dir, dir );
 690                         CrossProduct( normal, dir, hullDirs[j] );
 691                 }
 692
 693                 outside = qfalse;
 694                 for ( j = 0 ; j < numHullPoints ; j++ ) {
 695                         VectorSubtract( p, hullPoints[j], dir );
 696                         d = DotProduct( dir, hullDirs[j] );
 697                         if ( d >= ON_EPSILON ) {
 698                                 outside = qtrue;
 699                         }
 700                         if ( d >= -ON_EPSILON ) {
 701                                 hullSide[j] = qtrue;
 702                         } else {
 703                                 hullSide[j] = qfalse;
 704                         }
 705                 }
 706
 707                 // if the point is effectively inside, do nothing
 708                 if ( !outside ) {
 709                         continue;
 710                 }
 711
 712                 // find the back side to front side transition
 713                 for ( j = 0 ; j < numHullPoints ; j++ ) {
 714                         if ( !hullSide[ j % numHullPoints ] && hullSide[ (j + 1) % numHullPoints ] ) {
 715                                 break;
 716                         }
 717                 }
 718                 if ( j == numHullPoints ) {
 719                         continue;
 720                 }
 721
 722                 // insert the point here
 723                 VectorCopy( p, newHullPoints[0] );
 724                 numNew = 1;
 725
 726                 // copy over all points that aren't double fronts
 727                 j = (j+1)%numHullPoints;
 728                 for ( k = 0 ; k < numHullPoints ; k++ ) {
 729                         if ( hullSide[ (j+k) % numHullPoints ] && hullSide[ (j+k+1) % numHullPoints ] ) {
 730                                 continue;
 731                         }
 732                         copy = hullPoints[ (j+k+1) % numHullPoints ];
 733                         VectorCopy( copy, newHullPoints[numNew] );
 734                         numNew++;
 735                 }
 736
 737                 numHullPoints = numNew;
 738                 memcpy( hullPoints, newHullPoints, numHullPoints * sizeof(vec3_t) );
 739         }
 740
 741         FreeWinding( *hull );
 742         w = AllocWinding( numHullPoints );
 743         w->numpoints = numHullPoints;
 744         *hull = w;
 745         memcpy( w->p, hullPoints, numHullPoints * sizeof(vec3_t) );
 746 }
 747
 748