gfx/povray/polygon.c

   1 /****************************************************************************
   2 *                   polygon.c
   3 *
   4 *  This module implements functions that manipulate polygons.
   5 *
   6 *  This module was written by Dieter Bayer [DB].
   7 *
   8 *  from Persistence of Vision(tm) Ray Tracer
   9 *  Copyright 1996,1999 Persistence of Vision Team
  10 *---------------------------------------------------------------------------
  11 *  NOTICE: This source code file is provided so that users may experiment
  12 *  with enhancements to POV-Ray and to port the software to platforms other
  13 *  than those supported by the POV-Ray Team.  There are strict rules under
  14 *  which you are permitted to use this file.  The rules are in the file
  15 *  named POVLEGAL.DOC which should be distributed with this file.
  16 *  If POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  17 *  Team Coordinator by email to team-coord@povray.org or visit us on the web at
  18 *  http://www.povray.org. The latest version of POV-Ray may be found at this site.
  19 *
  20 * This program is based on the popular DKB raytracer version 2.12.
  21 * DKBTrace was originally written by David K. Buck.
  22 * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  23 *
  24 *****************************************************************************/
  25
  26 /****************************************************************************
  27 *
  28 *  Explanation:
  29 *
  30 *  Syntax:
  31 *
  32 *  ---
  33 *
  34 *  The "inside polygon"-test was taken from:
  35 *
  36 *    E. Haines, "An Introduction to Ray-Tracing", ..., pp. 53-59
  37 *
  38 *  ---
  39 *
  40 *  May 1994 : Creation. [DB]
  41 *
  42 *  Oct 1994 : Changed polygon structure. Polygon points are now stored
  43 *             in a seperate structure. This - together with the use of
  44 *             a transformation - allows to keep just one point definition
  45 *             for multiple copies of one polygon. [DB]
  46 *
  47 *****************************************************************************/
  48
  49 #include "frame.h"
  50 #include "vector.h"
  51 #include "povproto.h"
  52 #include "matrices.h"
  53 #include "objects.h"
  54 #include "polygon.h"
  55 #include "povray.h"
  56
  57
  58
  59
  60 /*****************************************************************************
  61 * Local preprocessor defines
  62 ******************************************************************************/
  63
  64 /* Minimal intersection depth for a valid intersection. */
  65
  66 #define DEPTH_TOLERANCE 1.0e-8
  67
  68 /* If |x| < ZERO_TOLERANCE x is assumed to be 0. */
  69
  70 #define ZERO_TOLERANCE 1.0e-10
  71
  72
  73
  74
  75 /*****************************************************************************
  76 * Static functions
  77 ******************************************************************************/
  78
  79 static int intersect_poylgon (RAY *Ray, POLYGON *Polyg, DBL *Depth);
  80 static int in_polygon (int Number, UV_VECT *Points, DBL u, DBL v);
  81 static int  All_Polygon_Intersections (OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack);
  82 static int  Inside_Polygon (VECTOR point, OBJECT *Object);
  83 static void Polygon_Normal (VECTOR Result, OBJECT *Object, INTERSECTION *Inter);
  84 static POLYGON *Copy_Polygon (OBJECT *Object);
  85 static void Translate_Polygon (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
  86 static void Rotate_Polygon (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
  87 static void Scale_Polygon (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
  88 static void Transform_Polygon (OBJECT *Object, TRANSFORM *Trans);
  89 static void Invert_Polygon (OBJECT *Object);
  90 static void Destroy_Polygon (OBJECT *Object);
  91 static void Compute_Polygon_BBox (POLYGON *Polyg);
  92
  93
  94
  95 /*****************************************************************************
  96 * Local variables
  97 ******************************************************************************/
  98
  99 static METHODS Polygon_Methods =
 100 {
 101   All_Polygon_Intersections,
 102   Inside_Polygon, Polygon_Normal,
 103   (COPY_METHOD)Copy_Polygon,
 104   Translate_Polygon, Rotate_Polygon,
 105   Scale_Polygon, Transform_Polygon, Invert_Polygon, Destroy_Polygon
 106 };
 107
 108
 109
 110 /*****************************************************************************
 111 *
 112 * FUNCTION
 113 *
 114 *   All_Polygon_Intersections
 115 *
 116 * INPUT
 117 *
 118 *   Object      - Object
 119 *   Ray         - Ray
 120 *   Depth_Stack - Intersection stack
 121 *
 122 * OUTPUT
 123 *
 124 *   Depth_Stack
 125 *
 126 * RETURNS
 127 *
 128 *   int - TRUE, if a intersection was found
 129 *
 130 * AUTHOR
 131 *
 132 *   Dieter Bayer
 133 *
 134 * DESCRIPTION
 135 *
 136 *   Determine ray/polygon intersection and clip intersection found.
 137 *
 138 * CHANGES
 139 *
 140 *   May 1994 : Creation.
 141 *
 142 ******************************************************************************/
 143
 144 static int All_Polygon_Intersections(OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack)
 145 {
 146   DBL Depth;
 147   VECTOR IPoint;
 148
 149   if (intersect_poylgon(Ray, (POLYGON *)Object, &Depth))
 150   {
 151     VEvaluateRay(IPoint, Ray->Initial, Depth, Ray->Direction);
 152
 153     if (Point_In_Clip(IPoint, Object->Clip))
 154     {
 155       push_entry(Depth, IPoint, Object, Depth_Stack);
 156
 157       return(TRUE);
 158     }
 159   }
 160
 161   return(FALSE);
 162 }
 163
 164
 165
 166 /*****************************************************************************
 167 *
 168 * FUNCTION
 169 *
 170 *   intersect_poylgon
 171 *
 172 * INPUT
 173 *
 174 *   Ray     - Ray
 175 *   Polyg - Polygon
 176 *   Depth   - Depth of intersection found
 177 *
 178 * OUTPUT
 179 *
 180 *   Depth
 181 *
 182 * RETURNS
 183 *
 184 *   int - TRUE if intersection found
 185 *
 186 * AUTHOR
 187 *
 188 *   Dieter Bayer
 189 *
 190 * DESCRIPTION
 191 *
 192 *   Determine ray/polygon intersection.
 193 *
 194 * CHANGES
 195 *
 196 *   May 1994 : Creation.
 197 *
 198 ******************************************************************************/
 199
 200 static int intersect_poylgon(RAY *Ray, POLYGON *Polyg, DBL *Depth)
 201 {
 202   DBL x, y, len;
 203   VECTOR p, d;
 204
 205   /* Don't test degenerate polygons. */
 206
 207   if (Test_Flag(Polyg, DEGENERATE_FLAG))
 208   {
 209     return(FALSE);
 210   }
 211
 212   Increase_Counter(stats[Ray_Polygon_Tests]);
 213
 214   /* Transform the ray into the polygon space. */
 215
 216   MInvTransPoint(p, Ray->Initial, Polyg->Trans);
 217
 218   MInvTransDirection(d, Ray->Direction, Polyg->Trans);
 219
 220   VLength(len, d);
 221
 222   VInverseScaleEq(d, len);
 223
 224   /* Intersect ray with the plane in which the polygon lies. */
 225
 226   if (fabs(d[Z]) < ZERO_TOLERANCE)
 227   {
 228     return(FALSE);
 229   }
 230
 231   *Depth = -p[Z] / d[Z];
 232
 233   if ((*Depth < DEPTH_TOLERANCE) || (*Depth > Max_Distance))
 234   {
 235     return(FALSE);
 236   }
 237
 238   /* Does the intersection point lie inside the polygon? */
 239
 240   x = p[X] + *Depth * d[X];
 241   y = p[Y] + *Depth * d[Y];
 242
 243   if (in_polygon(Polyg->Data->Number, Polyg->Data->Points, x, y))
 244   {
 245     Increase_Counter(stats[Ray_Polygon_Tests_Succeeded]);
 246
 247     *Depth /= len;
 248
 249     return (TRUE);
 250   }
 251   else
 252   {
 253     return (FALSE);
 254   }
 255 }
 256
 257
 258
 259 /*****************************************************************************
 260 *
 261 * FUNCTION
 262 *
 263 *   Inside_Polygon
 264 *
 265 * INPUT
 266 *
 267 *   IPoint - Intersection point
 268 *   Object - Object
 269 *
 270 * OUTPUT
 271 *
 272 * RETURNS
 273 *
 274 *   int - always FALSE
 275 *
 276 * AUTHOR
 277 *
 278 *   Dieter Bayer
 279 *
 280 * DESCRIPTION
 281 *
 282 *   Polygons can't be used in CSG.
 283 *
 284 * CHANGES
 285 *
 286 *   May 1994 : Creation.
 287 *
 288 ******************************************************************************/
 289
 290 static int Inside_Polygon(VECTOR IPoint, OBJECT *Object)
 291 {
 292   return(FALSE);
 293 }
 294
 295
 296
 297 /*****************************************************************************
 298 *
 299 * FUNCTION
 300 *
 301 *   Polygon_Normal
 302 *
 303 * INPUT
 304 *
 305 *   Result - Normal vector
 306 *   Object - Object
 307 *   Inter  - Intersection found
 308 *
 309 * OUTPUT
 310 *
 311 *   Result
 312 *
 313 * RETURNS
 314 *
 315 * AUTHOR
 316 *
 317 *   Dieter Bayer
 318 *
 319 * DESCRIPTION
 320 *
 321 *   Calculate the normal of the polygon in a given point.
 322 *
 323 * CHANGES
 324 *
 325 *   May 1994 : Creation.
 326 *
 327 ******************************************************************************/
 328
 329 static void Polygon_Normal(VECTOR Result, OBJECT *Object, INTERSECTION *Inter)
 330 {
 331   Assign_Vector(Result, ((POLYGON *)Object)->S_Normal);
 332 }
 333
 334
 335
 336 /*****************************************************************************
 337 *
 338 * FUNCTION
 339 *
 340 *   Translate_Polygon
 341 *
 342 * INPUT
 343 *
 344 *   Object - Object
 345 *   Vector - Translation vector
 346 *
 347 * OUTPUT
 348 *
 349 *   Object
 350 *
 351 * RETURNS
 352 *
 353 * AUTHOR
 354 *
 355 *   Dieter Bayer
 356 *
 357 * DESCRIPTION
 358 *
 359 *   Translate a polygon.
 360 *
 361 * CHANGES
 362 *
 363 *   May 1994 : Creation.
 364 *
 365 ******************************************************************************/
 366
 367 static void Translate_Polygon(OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
 368 {
 369   Transform_Polygon(Object, Trans);
 370 }
 371
 372
 373
 374 /*****************************************************************************
 375 *
 376 * FUNCTION
 377 *
 378 *   Rotate_Polygon
 379 *
 380 * INPUT
 381 *
 382 *   Object - Object
 383 *   Vector - Rotation vector
 384 *
 385 * OUTPUT
 386 *
 387 *   Object
 388 *
 389 * RETURNS
 390 *
 391 * AUTHOR
 392 *
 393 *   Dieter Bayer
 394 *
 395 * DESCRIPTION
 396 *
 397 *   Rotate a polygon.
 398 *
 399 * CHANGES
 400 *
 401 *   May 1994 : Creation.
 402 *
 403 ******************************************************************************/
 404
 405 static void Rotate_Polygon(OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
 406 {
 407   Transform_Polygon(Object, Trans);
 408 }
 409
 410
 411
 412 /*****************************************************************************
 413 *
 414 * FUNCTION
 415 *
 416 *   Scale_Polygon
 417 *
 418 * INPUT
 419 *
 420 *   Object - Object
 421 *   Vector - Scaling vector
 422 *
 423 * OUTPUT
 424 *
 425 *   Object
 426 *
 427 * RETURNS
 428 *
 429 * AUTHOR
 430 *
 431 *   Dieter Bayer
 432 *
 433 * DESCRIPTION
 434 *
 435 *   Scale a polygon.
 436 *
 437 * CHANGES
 438 *
 439 *   May 1994 : Creation.
 440 *
 441 ******************************************************************************/
 442
 443 static void Scale_Polygon(OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
 444 {
 445   Transform_Polygon(Object, Trans);
 446 }
 447
 448
 449
 450 /*****************************************************************************
 451 *
 452 * FUNCTION
 453 *
 454 *   Transform_Polygon
 455 *
 456 * INPUT
 457 *
 458 *   Object - Object
 459 *   Trans  - Transformation to apply
 460 *
 461 * OUTPUT
 462 *
 463 *   Object
 464 *
 465 * RETURNS
 466 *
 467 * AUTHOR
 468 *
 469 *   Dieter Bayer
 470 *
 471 * DESCRIPTION
 472 *
 473 *   Transform a polygon by transforming all points
 474 *   and recalculating the polygon.
 475 *
 476 * CHANGES
 477 *
 478 *   May 1994 : Creation.
 479 *
 480 ******************************************************************************/
 481
 482 static void Transform_Polygon(OBJECT *Object, TRANSFORM *Trans)
 483 {
 484   VECTOR N;
 485   POLYGON *Polyg = (POLYGON *)Object;
 486
 487   Compose_Transforms(Polyg->Trans, Trans);
 488
 489   Make_Vector(N, 0.0, 0.0, 1.0);
 490   MTransNormal(Polyg->S_Normal, N, Polyg->Trans);
 491
 492   VNormalizeEq(Polyg->S_Normal);
 493
 494   Compute_Polygon_BBox(Polyg);
 495 }
 496
 497
 498
 499 /*****************************************************************************
 500 *
 501 * FUNCTION
 502 *
 503 *   Invert_Polygon
 504 *
 505 * INPUT
 506 *
 507 *   Object - Object
 508 *
 509 * OUTPUT
 510 *
 511 *   Object
 512 *
 513 * RETURNS
 514 *
 515 * AUTHOR
 516 *
 517 *   Dieter Bayer
 518 *
 519 * DESCRIPTION
 520 *
 521 *   Invert a polygon.
 522 *
 523 * CHANGES
 524 *
 525 *   May 1994 : Creation.
 526 *
 527 ******************************************************************************/
 528
 529 static void Invert_Polygon(OBJECT *Object)
 530 {
 531 }
 532
 533
 534
 535 /*****************************************************************************
 536 *
 537 * FUNCTION
 538 *
 539 *   Create_Polygon
 540 *
 541 * INPUT
 542 *
 543 * OUTPUT
 544 *
 545 * RETURNS
 546 *
 547 *   POLYGON * - new polygon
 548 *
 549 * AUTHOR
 550 *
 551 *   Dieter Bayer
 552 *
 553 * DESCRIPTION
 554 *
 555 *   Create a new polygon.
 556 *
 557 * CHANGES
 558 *
 559 *   May 1994 : Creation.
 560 *
 561 ******************************************************************************/
 562
 563 POLYGON *Create_Polygon()
 564 {
 565   POLYGON *New;
 566
 567   New = (POLYGON *)POV_MALLOC(sizeof(POLYGON), "polygon");
 568
 569   INIT_OBJECT_FIELDS(New,POLYGON_OBJECT,&Polygon_Methods)
 570
 571   New->Trans = Create_Transform();
 572
 573   Make_Vector(New->S_Normal, 0.0, 0.0, 1.0);
 574
 575   New->Data = NULL;
 576
 577   return (New);
 578 }
 579
 580
 581
 582 /*****************************************************************************
 583 *
 584 * FUNCTION
 585 *
 586 *   Copy_Polygon
 587 *
 588 * INPUT
 589 *
 590 *   Object - Object
 591 *
 592 * OUTPUT
 593 *
 594 * RETURNS
 595 *
 596 *   void * - New polygon
 597 *
 598 * AUTHOR
 599 *
 600 *   Dieter Bayer
 601 *
 602 * DESCRIPTION
 603 *
 604 *   Copy a polygon structure.
 605 *
 606 * CHANGES
 607 *
 608 *   May 1994 : Creation.
 609 *
 610 ******************************************************************************/
 611
 612 static POLYGON *Copy_Polygon(OBJECT *Object)
 613 {
 614   POLYGON *New, *Polyg = (POLYGON *)Object;
 615
 616   New = Create_Polygon ();
 617
 618   /* Get rid of the transformation created in Create_Polygon(). */
 619
 620   Destroy_Transform(New->Trans);
 621
 622   /* Copy polygon. */
 623
 624   *New = *Polyg;
 625
 626   New->Trans = Copy_Transform(Polyg->Trans);
 627
 628   New->Data->References++;
 629
 630   return (New);
 631 }
 632
 633
 634
 635 /*****************************************************************************
 636 *
 637 * FUNCTION
 638 *
 639 *   Destroy_Polygon
 640 *
 641 * INPUT
 642 *
 643 *   Object - Object
 644 *
 645 * OUTPUT
 646 *
 647 *   Object
 648 *
 649 * RETURNS
 650 *
 651 * AUTHOR
 652 *
 653 *   Dieter Bayer
 654 *
 655 * DESCRIPTION
 656 *
 657 *   Destroy a polygon.
 658 *
 659 * CHANGES
 660 *
 661 *   May 1994 : Creation.
 662 *
 663 *   Dec 1994 : Fixed memory leakage. [DB]
 664 *
 665 ******************************************************************************/
 666
 667 static void Destroy_Polygon(OBJECT *Object)
 668 {
 669   POLYGON *Polyg = (POLYGON *)Object;
 670
 671   if (--(Polyg->Data->References) == 0)
 672   {
 673     POV_FREE (Polyg->Data->Points);
 674
 675     POV_FREE (Polyg->Data);
 676   }
 677
 678   Destroy_Transform(Polyg->Trans);
 679
 680   POV_FREE (Object);
 681 }
 682
 683
 684
 685 /*****************************************************************************
 686 *
 687 * FUNCTION
 688 *
 689 *   Compute_Polygon
 690 *
 691 * INPUT
 692 *
 693 *   Polyg - Polygon
 694 *   Points  - 3D points describing the polygon
 695 *
 696 * OUTPUT
 697 *
 698 *   Polyg
 699 *
 700 * RETURNS
 701 *
 702 * AUTHOR
 703 *
 704 *   Dieter Bayer
 705 *
 706 * DESCRIPTION
 707 *
 708 *   Compute the following things for a given polygon:
 709 *
 710 *     - Polygon transformation
 711 *
 712 *     - Array of 2d points describing the shape of the polygon
 713 *
 714 * CHANGES
 715 *
 716 *   May 1994 : Creation.
 717 *
 718 ******************************************************************************/
 719
 720 void Compute_Polygon(POLYGON *Polyg, int Number, VECTOR *Points)
 721 {
 722   int i;
 723   DBL x, y, z, d;
 724   VECTOR o, u = {}, v = {}, w = {}, N;
 725   MATRIX a, b;
 726
 727   /* Create polygon data. */
 728
 729   if (Polyg->Data == NULL)
 730   {
 731     Polyg->Data = (POLYGON_DATA *)POV_MALLOC(sizeof(POLYGON_DATA), "polygon points");
 732
 733     Polyg->Data->References = 1;
 734
 735     Polyg->Data->Number = Number;
 736
 737     Polyg->Data->Points = (UV_VECT *)POV_MALLOC(Number*sizeof(UV_VECT), "polygon points");
 738   }
 739   else
 740   {
 741     Error("Polygon data already computed.");
 742   }
 743
 744   /* Get polygon's coordinate system (one of the many possible) */
 745
 746   Assign_Vector(o, Points[0]);
 747
 748   /* Find valid, i.e. non-zero u vector. */
 749
 750   for (i = 1; i < Number; i++)
 751   {
 752     VSub(u, Points[i], o);
 753
 754     if (VSumSqr(u) > EPSILON)
 755     {
 756       break;
 757     }
 758   }
 759
 760   if (i == Number)
 761   {
 762     Set_Flag(Polyg, DEGENERATE_FLAG);
 763
 764     Warning(0.0, "Points in polygon are co-linear. Ignoring polygon.\n");
 765   }
 766
 767   /* Find valid, i.e. non-zero v and w vectors. */
 768
 769   for (i++; i < Number; i++)
 770   {
 771     VSub(v, Points[i], o);
 772
 773     VCross(w, u, v);
 774
 775     if ((VSumSqr(v) > EPSILON) && (VSumSqr(w) > EPSILON))
 776     {
 777       break;
 778     }
 779   }
 780
 781   if (i == Number)
 782   {
 783     Set_Flag(Polyg, DEGENERATE_FLAG);
 784
 785     Warning(0.0, "Points in polygon are co-linear. Ignoring polygon.\n");
 786   }
 787
 788   VCross(u, v, w);
 789   VCross(v, w, u);
 790
 791   VNormalize(u, u);
 792   VNormalize(v, v);
 793   VNormalize(w, w);
 794
 795   MIdentity(a);
 796   MIdentity(b);
 797
 798   a[3][0] = -o[X];
 799   a[3][1] = -o[Y];
 800   a[3][2] = -o[Z];
 801
 802   b[0][0] =  u[X];
 803   b[1][0] =  u[Y];
 804   b[2][0] =  u[Z];
 805
 806   b[0][1] =  v[X];
 807   b[1][1] =  v[Y];
 808   b[2][1] =  v[Z];
 809
 810   b[0][2] =  w[X];
 811   b[1][2] =  w[Y];
 812   b[2][2] =  w[Z];
 813
 814   MTimes(Polyg->Trans->inverse, a, b);
 815
 816   MInvers(Polyg->Trans->matrix, Polyg->Trans->inverse);
 817
 818   /* Project points onto the u,v-plane (3D --> 2D) */
 819
 820   for (i = 0; i < Number; i++)
 821   {
 822     x = Points[i][X] - o[X];
 823     y = Points[i][Y] - o[Y];
 824     z = Points[i][Z] - o[Z];
 825
 826     d = x * w[X] + y * w[Y] + z * w[Z];
 827
 828     if (fabs(d) > ZERO_TOLERANCE)
 829     {
 830       Set_Flag(Polyg, DEGENERATE_FLAG);
 831
 832       Warning(0.0, "Points in polygon are not co-planar. Ignoring polygons.\n");
 833     }
 834
 835     Polyg->Data->Points[i][X] = x * u[X] + y * u[Y] + z * u[Z];
 836     Polyg->Data->Points[i][Y] = x * v[X] + y * v[Y] + z * v[Z];
 837   }
 838
 839   Make_Vector(N, 0.0, 0.0, 1.0);
 840   MTransNormal(Polyg->S_Normal, N, Polyg->Trans);
 841
 842   VNormalizeEq(Polyg->S_Normal);
 843
 844   Compute_Polygon_BBox(Polyg);
 845 }
 846
 847
 848
 849 /*****************************************************************************
 850 *
 851 * FUNCTION
 852 *
 853 *   Compute_Polygon_BBox
 854 *
 855 * INPUT
 856 *
 857 *   Polyg - Polygon
 858 *
 859 * OUTPUT
 860 *
 861 *   Polyg
 862 *
 863 * RETURNS
 864 *
 865 * AUTHOR
 866 *
 867 *   Dieter Bayer
 868 *
 869 * DESCRIPTION
 870 *
 871 *   Calculate the bounding box of a polygon.
 872 *
 873 * CHANGES
 874 *
 875 *   May 1994 : Creation.
 876 *
 877 ******************************************************************************/
 878
 879 static void Compute_Polygon_BBox(POLYGON *Polyg)
 880 {
 881   int i;
 882   VECTOR p, Puv, Min, Max;
 883
 884   Min[X] = Min[Y] = Min[Z] =  BOUND_HUGE;
 885   Max[X] = Max[Y] = Max[Z] = -BOUND_HUGE;
 886
 887   for (i = 0; i < Polyg->Data->Number; i++)
 888   {
 889     Puv[X] = Polyg->Data->Points[i][X];
 890     Puv[Y] = Polyg->Data->Points[i][Y];
 891     Puv[Z] = 0.0;
 892
 893     MTransPoint(p, Puv, Polyg->Trans);
 894
 895     Min[X] = min(Min[X], p[X]);
 896     Min[Y] = min(Min[Y], p[Y]);
 897     Min[Z] = min(Min[Z], p[Z]);
 898     Max[X] = max(Max[X], p[X]);
 899     Max[Y] = max(Max[Y], p[Y]);
 900     Max[Z] = max(Max[Z], p[Z]);
 901   }
 902
 903   Make_BBox_from_min_max(Polyg->BBox, Min, Max);
 904
 905   if (fabs(Polyg->BBox.Lengths[X]) < Small_Tolerance)
 906   {
 907     Polyg->BBox.Lower_Left[X] -= Small_Tolerance;
 908     Polyg->BBox.Lengths[X]    += 2.0 * Small_Tolerance;
 909   }
 910
 911   if (fabs(Polyg->BBox.Lengths[Y]) < Small_Tolerance)
 912   {
 913     Polyg->BBox.Lower_Left[Y] -= Small_Tolerance;
 914     Polyg->BBox.Lengths[Y]    += 2.0 * Small_Tolerance;
 915   }
 916
 917   if (fabs(Polyg->BBox.Lengths[Z]) < Small_Tolerance)
 918   {
 919     Polyg->BBox.Lower_Left[Z] -= Small_Tolerance;
 920     Polyg->BBox.Lengths[Z]    += 2.0 * Small_Tolerance;
 921   }
 922 }
 923
 924
 925
 926 /*****************************************************************************
 927 *
 928 * FUNCTION
 929 *
 930 *   in_polygon
 931 *
 932 * INPUT
 933 *
 934 *   Number - Number of points
 935 *   Points - Points describing polygon's shape
 936 *   u, v   - 2D-coordinates of the point to test
 937 *
 938 * OUTPUT
 939 *
 940 * RETURNS
 941 *
 942 *   int - TRUE, if inside
 943 *
 944 * AUTHOR
 945 *
 946 *   Eric Haines, 3D/Eye Inc, erich@eye.com
 947 *
 948 * DESCRIPTION
 949 *
 950 * ======= Crossings Multiply algorithm ===================================
 951 *
 952 * This version is usually somewhat faster than the original published in
 953 * Graphics Gems IV; by turning the division for testing the X axis crossing
 954 * into a tricky multiplication test this part of the test became faster,
 955 * which had the additional effect of making the test for "both to left or
 956 * both to right" a bit slower for triangles than simply computing the
 957 * intersection each time.  The main increase is in triangle testing speed,
 958 * which was about 15% faster; all other polygon complexities were pretty much
 959 * the same as before.  On machines where division is very expensive (not the
 960 * case on the HP 9000 series on which I tested) this test should be much
 961 * faster overall than the old code.  Your mileage may (in fact, will) vary,
 962 * depending on the machine and the test data, but in general I believe this
 963 * code is both shorter and faster.  This test was inspired by unpublished
 964 * Graphics Gems submitted by Joseph Samosky and Mark Haigh-Hutchinson.
 965 * Related work by Samosky is in:
 966 *
 967 * Samosky, Joseph, "SectionView: A system for interactively specifying and
 968 * visualizing sections through three-dimensional medical image data",
 969 * M.S. Thesis, Department of Electrical Engineering and Computer Science,
 970 * Massachusetts Institute of Technology, 1993.
 971 *
 972 *
 973 * Shoot a test ray along +X axis.  The strategy is to compare vertex Y values
 974 * to the testing point's Y and quickly discard edges which are entirely to one
 975 * side of the test ray.
 976 *
 977 * CHANGES
 978 *
 979 *   -
 980 *
 981 ******************************************************************************/
 982
 983 static int in_polygon(int Number, UV_VECT *Points, DBL u, DBL  v)
 984 {
 985   register int i, yflag0, yflag1, inside_flag;
 986   register DBL ty, tx;
 987   register DBL *vtx0, *vtx1, *first;
 988
 989   tx = u;
 990   ty = v;
 991
 992   vtx0 = &Points[0][X];
 993   vtx1 = &Points[1][X];
 994
 995   first = vtx0;
 996
 997   /* get test bit for above/below X axis */
 998
 999   yflag0 = (vtx0[Y] >= ty);
1000
1001   inside_flag = FALSE;
1002
1003   for (i = 1; i < Number; )
1004   {
1005     yflag1 = (vtx1[Y] >= ty);
1006
1007     /*
1008      * Check if endpoints straddle (are on opposite sides) of X axis
1009      * (i.e. the Y's differ); if so, +X ray could intersect this edge.
1010      * The old test also checked whether the endpoints are both to the
1011      * right or to the left of the test point.  However, given the faster
1012      * intersection point computation used below, this test was found to
1013      * be a break-even proposition for most polygons and a loser for
1014      * triangles (where 50% or more of the edges which survive this test
1015      * will cross quadrants and so have to have the X intersection computed
1016      * anyway).  I credit Joseph Samosky with inspiring me to try dropping
1017      * the "both left or both right" part of my code.
1018      */
1019
1020     if (yflag0 != yflag1)
1021     {
1022       /*
1023        * Check intersection of pgon segment with +X ray.
1024        * Note if >= point's X; if so, the ray hits it.
1025        * The division operation is avoided for the ">=" test by checking
1026        * the sign of the first vertex wrto the test point; idea inspired
1027        * by Joseph Samosky's and Mark Haigh-Hutchinson's different
1028        * polygon inclusion tests.
1029        */
1030
1031       if (((vtx1[Y]-ty) * (vtx0[X]-vtx1[X]) >= (vtx1[X]-tx) * (vtx0[Y]-vtx1[Y])) == yflag1)
1032       {
1033         inside_flag = !inside_flag;
1034       }
1035     }
1036
1037     /* Move to the next pair of vertices, retaining info as possible. */
1038
1039     if ((i < Number-2) && (vtx1[X] == first[X]) && (vtx1[Y] == first[Y]))
1040     {
1041       vtx0 = &Points[++i][X];
1042       vtx1 = &Points[++i][X];
1043
1044       yflag0 = (vtx0[Y] >= ty);
1045
1046       first = vtx0;
1047     }
1048     else
1049     {
1050       vtx0 = vtx1;
1051       vtx1 = &Points[++i][X];
1052
1053       yflag0 = yflag1;
1054     }
1055   }
1056
1057   return(inside_flag);
1058 }
1059