tools/quake2/qdata_heretic2/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 "inout.h"
  25 #include "mathlib.h"
  26 #include "polylib.h"
  27
  28
  29 extern int numthreads;
  30
  31 // counters are only bumped when running single threaded,
  32 // because they are an awefull coherence problem
  33 int c_active_windings;
  34 int c_peak_windings;
  35 int c_winding_allocs;
  36 int c_winding_points;
  37
  38 #define BOGUS_RANGE 8192
  39
  40 void pw( winding_t *w ){
  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         winding_t   *w;
  54         int s;
  55
  56         if ( numthreads == 1 ) {
  57                 c_winding_allocs++;
  58                 c_winding_points += points;
  59                 c_active_windings++;
  60                 if ( c_active_windings > c_peak_windings ) {
  61                         c_peak_windings = c_active_windings;
  62                 }
  63         }
  64         s = sizeof( vec_t ) * 3 * points + sizeof( int );
  65         w = malloc( s );
  66         if ( !w ) {
  67                 Error( "AllocWinding MALLOC failed!  Could not allocate %s bytes.", s );
  68         }
  69         memset( w, 0, s );
  70         return w;
  71 }
  72
  73 void FreeWinding( winding_t *w ){
  74         if ( *(unsigned *)w == 0xdeaddead ) {
  75                 Error( "FreeWinding: freed a freed winding" );
  76         }
  77         *(unsigned *)w = 0xdeaddead;
  78
  79         if ( numthreads == 1 ) {
  80                 c_active_windings--;
  81         }
  82         free( w );
  83 }
  84
  85 /*
  86    ============
  87    RemoveColinearPoints
  88    ============
  89  */
  90 int c_removed;
  91
  92 void    RemoveColinearPoints( winding_t *w ){
  93         int i, j, k;
  94         vec3_t v1, v2;
  95         int nump;
  96         vec3_t p[MAX_POINTS_ON_WINDING];
  97
  98         nump = 0;
  99         for ( i = 0 ; i < w->numpoints ; i++ )
 100         {
 101                 j = ( i + 1 ) % w->numpoints;
 102                 k = ( i + w->numpoints - 1 ) % w->numpoints;
 103                 VectorSubtract( w->p[j], w->p[i], v1 );
 104                 VectorSubtract( w->p[i], w->p[k], v2 );
 105                 VectorNormalize( v1,v1 );
 106                 VectorNormalize( v2,v2 );
 107                 if ( DotProduct( v1, v2 ) < 0.999 ) {
 108                         VectorCopy( w->p[i], p[nump] );
 109                         nump++;
 110                 }
 111         }
 112
 113         if ( nump == w->numpoints ) {
 114                 return;
 115         }
 116
 117         if ( numthreads == 1 ) {
 118                 c_removed += w->numpoints - nump;
 119         }
 120         w->numpoints = nump;
 121         memcpy( w->p, p, nump * sizeof( p[0] ) );
 122 }
 123
 124 /*
 125    ============
 126    WindingPlane
 127    ============
 128  */
 129 void WindingPlane( winding_t *w, vec3_t normal, vec_t *dist ){
 130         vec3_t v1, v2;
 131
 132         VectorSubtract( w->p[1], w->p[0], v1 );
 133         VectorSubtract( w->p[2], w->p[0], v2 );
 134         CrossProduct( v2, v1, normal );
 135         VectorNormalize( normal, normal );
 136         *dist = DotProduct( w->p[0], normal );
 137
 138 }
 139
 140 /*
 141    =============
 142    WindingArea
 143    =============
 144  */
 145 vec_t   WindingArea( winding_t *w ){
 146         int i;
 147         vec3_t d1, d2, cross;
 148         vec_t total;
 149
 150         total = 0;
 151         for ( i = 2 ; i < w->numpoints ; i++ )
 152         {
 153                 VectorSubtract( w->p[i - 1], w->p[0], d1 );
 154                 VectorSubtract( w->p[i], w->p[0], d2 );
 155                 CrossProduct( d1, d2, cross );
 156                 total += 0.5 * VectorLength( cross );
 157         }
 158         return total;
 159 }
 160
 161 void    WindingBounds( winding_t *w, vec3_t mins, vec3_t maxs ){
 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                         }
 176                         if ( v > maxs[j] ) {
 177                                 maxs[j] = v;
 178                         }
 179                 }
 180         }
 181 }
 182
 183 /*
 184    =============
 185    WindingCenter
 186    =============
 187  */
 188 void    WindingCenter( winding_t *w, vec3_t center ){
 189         int i;
 190         float scale;
 191
 192         VectorCopy( vec3_origin, center );
 193         for ( i = 0 ; i < w->numpoints ; i++ )
 194                 VectorAdd( w->p[i], center, center );
 195
 196         scale = 1.0 / w->numpoints;
 197         VectorScale( center, scale, center );
 198 }
 199
 200 /*
 201    =================
 202    BaseWindingForPlane
 203    =================
 204  */
 205 winding_t *BaseWindingForPlane( vec3_t normal, vec_t dist ){
 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                         x = i;
 220                         max = v;
 221                 }
 222         }
 223         if ( x == -1 ) {
 224                 Error( "BaseWindingForPlane: no axis found" );
 225         }
 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         int size;
 277         winding_t   *c;
 278
 279         c = AllocWinding( w->numpoints );
 280         size = (int)( (winding_t *)0 )->p[w->numpoints];
 281         memcpy( c, w, size );
 282         return c;
 283 }
 284
 285 /*
 286    ==================
 287    ReverseWinding
 288    ==================
 289  */
 290 winding_t   *ReverseWinding( winding_t *w ){
 291         int i;
 292         winding_t   *c;
 293
 294         c = AllocWinding( w->numpoints );
 295         for ( i = 0 ; i < w->numpoints ; i++ )
 296         {
 297                 VectorCopy( w->p[w->numpoints - 1 - i], c->p[i] );
 298         }
 299         c->numpoints = w->numpoints;
 300         return c;
 301 }
 302
 303
 304 /*
 305    =============
 306    ClipWindingEpsilon
 307    =============
 308  */
 309 void    ClipWindingEpsilon( winding_t *in, vec3_t normal, vec_t dist,
 310                                                         vec_t epsilon, winding_t **front, winding_t **back ){
 311         vec_t dists[MAX_POINTS_ON_WINDING + 4];
 312         int sides[MAX_POINTS_ON_WINDING + 4];
 313         int counts[3];
 314         vec_t dot;          // VC 4.2 optimizer bug if not static
 315         int i, j;
 316         vec_t   *p1, *p2;
 317         vec3_t mid;
 318         winding_t   *f, *b;
 319         int maxpts;
 320
 321         if ( in->numpoints >= MAX_POINTS_ON_WINDING - 4 ) {
 322                 Error( "ClipWinding: MAX_POINTS_ON_WINDING" );
 323         }
 324
 325         counts[0] = counts[1] = counts[2] = 0;
 326
 327 // determine sides for each point
 328         for ( i = 0 ; i < in->numpoints ; i++ )
 329         {
 330                 dot = DotProduct( in->p[i], normal );
 331                 dot -= dist;
 332                 dists[i] = dot;
 333                 if ( dot > epsilon ) {
 334                         sides[i] = SIDE_FRONT;
 335                 }
 336                 else if ( dot < -epsilon ) {
 337                         sides[i] = SIDE_BACK;
 338                 }
 339                 else
 340                 {
 341                         sides[i] = SIDE_ON;
 342                 }
 343                 counts[sides[i]]++;
 344         }
 345         sides[i] = sides[0];
 346         dists[i] = dists[0];
 347
 348         *front = *back = NULL;
 349
 350         if ( !counts[0] ) {
 351                 *back = CopyWinding( in );
 352                 return;
 353         }
 354         if ( !counts[1] ) {
 355                 *front = CopyWinding( in );
 356                 return;
 357         }
 358
 359         maxpts = in->numpoints + 4;   // cant use counts[0]+2 because
 360                                       // of fp grouping errors
 361
 362         *front = f = AllocWinding( maxpts );
 363         *back = b = AllocWinding( maxpts );
 364
 365         for ( i = 0 ; i < in->numpoints ; i++ )
 366         {
 367                 p1 = in->p[i];
 368
 369                 if ( sides[i] == SIDE_ON ) {
 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                         VectorCopy( p1, f->p[f->numpoints] );
 379                         f->numpoints++;
 380                 }
 381                 if ( sides[i] == SIDE_BACK ) {
 382                         VectorCopy( p1, b->p[b->numpoints] );
 383                         b->numpoints++;
 384                 }
 385
 386                 if ( sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i] ) {
 387                         continue;
 388                 }
 389
 390                 // generate a split point
 391                 p2 = in->p[( i + 1 ) % in->numpoints];
 392
 393                 dot = dists[i] / ( dists[i] - dists[i + 1] );
 394                 for ( j = 0 ; j < 3 ; j++ )
 395                 {   // avoid round off error when possible
 396                         if ( normal[j] == 1 ) {
 397                                 mid[j] = dist;
 398                         }
 399                         else if ( normal[j] == -1 ) {
 400                                 mid[j] = -dist;
 401                         }
 402                         else{
 403                                 mid[j] = p1[j] + dot * ( p2[j] - p1[j] );
 404                         }
 405                 }
 406
 407                 VectorCopy( mid, f->p[f->numpoints] );
 408                 f->numpoints++;
 409                 VectorCopy( mid, b->p[b->numpoints] );
 410                 b->numpoints++;
 411         }
 412
 413         if ( f->numpoints > maxpts || b->numpoints > maxpts ) {
 414                 Error( "ClipWinding: points exceeded estimate" );
 415         }
 416         if ( f->numpoints > MAX_POINTS_ON_WINDING || b->numpoints > MAX_POINTS_ON_WINDING ) {
 417                 Error( "ClipWinding: MAX_POINTS_ON_WINDING" );
 418         }
 419 }
 420
 421
 422 /*
 423    =============
 424    ChopWindingInPlace
 425    =============
 426  */
 427 void ChopWindingInPlace( winding_t **inout, vec3_t normal, vec_t dist, vec_t epsilon ){
 428         winding_t   *in;
 429         vec_t dists[MAX_POINTS_ON_WINDING + 4];
 430         int sides[MAX_POINTS_ON_WINDING + 4];
 431         int counts[3];
 432         vec_t dot;          // VC 4.2 optimizer bug if not static
 433         int i, j;
 434         vec_t   *p1, *p2;
 435         vec3_t mid;
 436         winding_t   *f;
 437         int maxpts;
 438
 439         in = *inout;
 440         counts[0] = counts[1] = counts[2] = 0;
 441
 442         if ( !in ) {
 443                 printf( "Warning: NULL passed to ChopWindingInPlace\n" );
 444                 return;
 445         }
 446         if ( in->numpoints >= MAX_POINTS_ON_WINDING - 4 ) {
 447                 Error( "ChopWinding: MAX_POINTS_ON_WINDING" );
 448         }
 449
 450 // determine sides for each point
 451         for ( i = 0 ; i < in->numpoints ; i++ )
 452         {
 453                 dot = DotProduct( in->p[i], normal );
 454                 dot -= dist;
 455                 dists[i] = dot;
 456                 if ( dot > epsilon ) {
 457                         sides[i] = SIDE_FRONT;
 458                 }
 459                 else if ( dot < -epsilon ) {
 460                         sides[i] = SIDE_BACK;
 461                 }
 462                 else
 463                 {
 464                         sides[i] = SIDE_ON;
 465                 }
 466                 counts[sides[i]]++;
 467         }
 468         sides[i] = sides[0];
 469         dists[i] = dists[0];
 470
 471         if ( !counts[0] ) {
 472                 FreeWinding( in );
 473                 *inout = NULL;
 474                 return;
 475         }
 476         if ( !counts[1] ) {
 477                 return;     // inout stays the same
 478
 479         }
 480         maxpts = in->numpoints + 4;   // cant use counts[0]+2 because
 481                                       // of fp grouping errors
 482
 483         f = AllocWinding( maxpts );
 484
 485         for ( i = 0 ; i < in->numpoints ; i++ )
 486         {
 487                 p1 = in->p[i];
 488
 489                 if ( sides[i] == SIDE_ON ) {
 490                         VectorCopy( p1, f->p[f->numpoints] );
 491                         f->numpoints++;
 492                         continue;
 493                 }
 494
 495                 if ( sides[i] == SIDE_FRONT ) {
 496                         VectorCopy( p1, f->p[f->numpoints] );
 497                         f->numpoints++;
 498                 }
 499
 500                 if ( sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i] ) {
 501                         continue;
 502                 }
 503
 504                 // generate a split point
 505                 p2 = in->p[( i + 1 ) % in->numpoints];
 506
 507                 dot = dists[i] / ( dists[i] - dists[i + 1] );
 508                 for ( j = 0 ; j < 3 ; j++ )
 509                 {   // avoid round off error when possible
 510                         if ( normal[j] == 1 ) {
 511                                 mid[j] = dist;
 512                         }
 513                         else if ( normal[j] == -1 ) {
 514                                 mid[j] = -dist;
 515                         }
 516                         else{
 517                                 mid[j] = p1[j] + dot * ( p2[j] - p1[j] );
 518                         }
 519                 }
 520
 521                 VectorCopy( mid, f->p[f->numpoints] );
 522                 f->numpoints++;
 523         }
 524
 525         if ( f->numpoints > maxpts ) {
 526                 Error( "ClipWinding: points exceeded estimate" );
 527         }
 528         if ( f->numpoints > MAX_POINTS_ON_WINDING ) {
 529                 Error( "ClipWinding: MAX_POINTS_ON_WINDING" );
 530         }
 531
 532         FreeWinding( in );
 533         *inout = f;
 534 }
 535
 536
 537 /*
 538    =================
 539    ChopWinding
 540
 541    Returns the fragment of in that is on the front side
 542    of the cliping plane.  The original is freed.
 543    =================
 544  */
 545 winding_t   *ChopWinding( winding_t *in, vec3_t normal, vec_t dist ){
 546         winding_t   *f, *b;
 547
 548         ClipWindingEpsilon( in, normal, dist, ON_EPSILON, &f, &b );
 549         FreeWinding( in );
 550         if ( b ) {
 551                 FreeWinding( b );
 552         }
 553         return f;
 554 }
 555
 556
 557 /*
 558    =================
 559    CheckWinding
 560
 561    =================
 562  */
 563 void CheckWinding( winding_t *w ){
 564         int i, j;
 565         vec_t   *p1, *p2;
 566         vec_t d, edgedist;
 567         vec3_t dir, edgenormal, facenormal;
 568         vec_t area;
 569         vec_t facedist;
 570
 571         if ( w->numpoints < 3 ) {
 572                 Error( "CheckWinding: %i points",w->numpoints );
 573         }
 574
 575         area = WindingArea( w );
 576         if ( area < 1 ) {
 577                 Error( "CheckWinding: %f area", area );
 578         }
 579
 580         WindingPlane( w, facenormal, &facedist );
 581
 582         for ( i = 0 ; i < w->numpoints ; i++ )
 583         {
 584                 p1 = w->p[i];
 585
 586                 for ( j = 0 ; j < 3 ; j++ )
 587                         if ( p1[j] > BOGUS_RANGE || p1[j] < -BOGUS_RANGE ) {
 588                                 Error( "CheckFace: BUGUS_RANGE: %f",p1[j] );
 589                         }
 590
 591                 j = i + 1 == w->numpoints ? 0 : i + 1;
 592
 593                 // check the point is on the face plane
 594                 d = DotProduct( p1, facenormal ) - facedist;
 595                 if ( d < -ON_EPSILON || d > ON_EPSILON ) {
 596                         Error( "CheckWinding: point off plane" );
 597                 }
 598
 599                 // check the edge isnt degenerate
 600                 p2 = w->p[j];
 601                 VectorSubtract( p2, p1, dir );
 602
 603                 if ( VectorLength( dir ) < ON_EPSILON ) {
 604                         Error( "CheckWinding: degenerate edge" );
 605                 }
 606
 607                 CrossProduct( facenormal, dir, edgenormal );
 608                 VectorNormalize( edgenormal, edgenormal );
 609                 edgedist = DotProduct( p1, edgenormal );
 610                 edgedist += ON_EPSILON;
 611
 612                 // all other points must be on front side
 613                 for ( j = 0 ; j < w->numpoints ; j++ )
 614                 {
 615                         if ( j == i ) {
 616                                 continue;
 617                         }
 618                         d = DotProduct( w->p[j], edgenormal );
 619                         if ( d > edgedist ) {
 620                                 Error( "CheckWinding: non-convex" );
 621                         }
 622                 }
 623         }
 624 }
 625
 626
 627 /*
 628    ============
 629    WindingOnPlaneSide
 630    ============
 631  */
 632 int     WindingOnPlaneSide( winding_t *w, vec3_t normal, vec_t dist ){
 633         qboolean front, back;
 634         int i;
 635         vec_t d;
 636
 637         front = false;
 638         back = false;
 639         for ( i = 0 ; i < w->numpoints ; i++ )
 640         {
 641                 d = DotProduct( w->p[i], normal ) - dist;
 642                 if ( d < -ON_EPSILON ) {
 643                         if ( front ) {
 644                                 return SIDE_CROSS;
 645                         }
 646                         back = true;
 647                         continue;
 648                 }
 649                 if ( d > ON_EPSILON ) {
 650                         if ( back ) {
 651                                 return SIDE_CROSS;
 652                         }
 653                         front = true;
 654                         continue;
 655                 }
 656         }
 657
 658         if ( back ) {
 659                 return SIDE_BACK;
 660         }
 661         if ( front ) {
 662                 return SIDE_FRONT;
 663         }
 664         return SIDE_ON;
 665 }