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