gfx/povray/super.c

   1 /****************************************************************************
   2 *                   super.c
   3 *
   4 *  This module implements functions that manipulate superellipsoids.
   5 *
   6 *  Original code written by Alexander Enzmann.
   7 *  Adaption to POV-Ray by Dieter Bayer [DB].
   8 *
   9 *  from Persistence of Vision(tm) Ray Tracer
  10 *  Copyright 1996,1999 Persistence of Vision Team
  11 *---------------------------------------------------------------------------
  12 *  NOTICE: This source code file is provided so that users may experiment
  13 *  with enhancements to POV-Ray and to port the software to platforms other
  14 *  than those supported by the POV-Ray Team.  There are strict rules under
  15 *  which you are permitted to use this file.  The rules are in the file
  16 *  named POVLEGAL.DOC which should be distributed with this file.
  17 *  If POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  18 *  Team Coordinator by email to team-coord@povray.org or visit us on the web at
  19 *  http://www.povray.org. The latest version of POV-Ray may be found at this site.
  20 *
  21 * This program is based on the popular DKB raytracer version 2.12.
  22 * DKBTrace was originally written by David K. Buck.
  23 * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  24 *
  25 *****************************************************************************/
  26
  27 /****************************************************************************
  28 *
  29 *  Explanation:
  30 *
  31 *    Superellipsoids are defined by the implicit equation
  32 *
  33 *      f(x,y,z) = (|x|^(2/e) + |y|^(2/e))^(e/n) + |z|^(2/n) - 1
  34 *
  35 *    Where e is the east/west exponent and n is the north/south exponent.
  36 *
  37 *    NOTE: The exponents are precomputed and stored in the Power element.
  38 *
  39 *    NOTE: Values of e and n that are close to degenerate (e.g.,
  40 *          below around 0.1) appear to give the root solver fits.
  41 *          Not sure quite where the problem lays just yet.
  42 *
  43 *  Syntax:
  44 *
  45 *    superellipsoid { <e, n> }
  46 *
  47 *
  48 *  ---
  49 *
  50 *  Oct 1994 : Creation.
  51 *
  52 *****************************************************************************/
  53
  54 #include "frame.h"
  55 #include "povray.h"
  56 #include "vector.h"
  57 #include "povproto.h"
  58 #include "bbox.h"
  59 #include "matrices.h"
  60 #include "objects.h"
  61 #include "super.h"
  62
  63
  64
  65 /*****************************************************************************
  66 * Local preprocessor defines
  67 ******************************************************************************/
  68
  69 /* Minimal intersection depth for a valid intersection. */
  70
  71 #define DEPTH_TOLERANCE 1.0e-4
  72
  73 /* If |x| < ZERO_TOLERANCE it is regarded to be 0. */
  74
  75 #define ZERO_TOLERANCE 1.0e-10
  76
  77 /* This is not the signum function because SGNX(0) is 1. */
  78
  79 #define SGNX(x) (((x) >= 0.0) ? 1.0 : -1.0)
  80
  81 #define MIN_VALUE -1.0
  82 #define MAX_VALUE  1.0
  83
  84 #define MAX_ITERATIONS 20
  85
  86 #define PLANECOUNT 9
  87
  88
  89
  90 /*****************************************************************************
  91 * Local typedefs
  92 ******************************************************************************/
  93
  94
  95
  96 /*****************************************************************************
  97 * Static functions
  98 ******************************************************************************/
  99
 100 static int intersect_superellipsoid (RAY *Ray, SUPERELLIPSOID *Superellipsoid, ISTACK *Depth_Stack);
 101 static int intersect_box (VECTOR P, VECTOR D, DBL *dmin, DBL *dmax);
 102 static DBL power (DBL x, DBL e);
 103 static DBL evaluate_superellipsoid (VECTOR P, SUPERELLIPSOID *Superellipsoid);
 104 static int compdists (CONST void *in_a, CONST void *in_b);
 105
 106 static int find_ray_plane_points (VECTOR P, VECTOR D, int cnt, DBL *dists, DBL mindist, DBL maxdist);
 107 static void solve_hit1 (SUPERELLIPSOID *Super, DBL v0, VECTOR tP0, DBL v1, VECTOR tP1, VECTOR P);
 108 static int check_hit2 (SUPERELLIPSOID *Super, VECTOR P, VECTOR D, DBL t0, VECTOR P0, DBL v0, DBL t1, DBL *t, VECTOR Q);
 109 static int insert_hit (OBJECT *Object, RAY *Ray, DBL Depth, ISTACK *Depth_Stack);
 110
 111 static int   All_Superellipsoid_Intersections (OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack);
 112 static int   Inside_Superellipsoid (VECTOR point, OBJECT *Object);
 113 static void  Superellipsoid_Normal (VECTOR Result, OBJECT *Object, INTERSECTION *Inter);
 114 static SUPERELLIPSOID *Copy_Superellipsoid (OBJECT *Object);
 115 static void  Translate_Superellipsoid (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
 116 static void  Rotate_Superellipsoid (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
 117 static void  Scale_Superellipsoid (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
 118 static void  Transform_Superellipsoid (OBJECT *Object, TRANSFORM *Trans);
 119 static void  Invert_Superellipsoid (OBJECT *Object);
 120 static void  Destroy_Superellipsoid (OBJECT *Object);
 121
 122 /*****************************************************************************
 123 * Local variables
 124 ******************************************************************************/
 125
 126 static METHODS Superellipsoid_Methods =
 127 {
 128   All_Superellipsoid_Intersections,
 129   Inside_Superellipsoid, Superellipsoid_Normal,
 130   (COPY_METHOD)Copy_Superellipsoid,
 131   Translate_Superellipsoid, Rotate_Superellipsoid,
 132   Scale_Superellipsoid, Transform_Superellipsoid, Invert_Superellipsoid, Destroy_Superellipsoid
 133 };
 134
 135 static DBL planes[PLANECOUNT][4] =
 136   {{1, 1, 0, 0}, {1,-1, 0, 0},
 137    {1, 0, 1, 0}, {1, 0,-1, 0},
 138    {0, 1, 1, 0}, {0, 1,-1, 0},
 139    {1, 0, 0, 0},
 140    {0, 1, 0, 0},
 141    {0, 0, 1, 0}};
 142
 143
 144
 145 /*****************************************************************************
 146 *
 147 * FUNCTION
 148 *
 149 *   All_Superellipsoid_Intersections
 150 *
 151 * INPUT
 152 *
 153 *   Object      - Object
 154 *   Ray         - Ray
 155 *   Depth_Stack - Intersection stack
 156 *
 157 * OUTPUT
 158 *
 159 *   Depth_Stack
 160 *
 161 * RETURNS
 162 *
 163 *   int - TRUE, if a intersection was found
 164 *
 165 * AUTHOR
 166 *
 167 *   Dieter Bayer
 168 *
 169 * DESCRIPTION
 170 *
 171 *   Determine ray/superellipsoid intersection and clip intersection found.
 172 *
 173 * CHANGES
 174 *
 175 *   Oct 1994 : Creation.
 176 *
 177 ******************************************************************************/
 178
 179 static int All_Superellipsoid_Intersections(OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack)
 180 {
 181   Increase_Counter(stats[Ray_Superellipsoid_Tests]);
 182
 183   if (intersect_superellipsoid(Ray, (SUPERELLIPSOID *)Object, Depth_Stack))
 184   {
 185     Increase_Counter(stats[Ray_Superellipsoid_Tests_Succeeded]);
 186
 187     return(TRUE);
 188   }
 189   else
 190   {
 191     return(FALSE);
 192   }
 193 }
 194
 195
 196
 197 /*****************************************************************************
 198 *
 199 * FUNCTION
 200 *
 201 *   intersect_superellipsoid
 202 *
 203 * INPUT
 204 *
 205 *   Ray            - Ray
 206 *   Superellipsoid - Superellipsoid
 207 *   Depth_Stack    - Depth stack
 208 *
 209 * OUTPUT
 210 *
 211 *   Intersection
 212 *
 213 * RETURNS
 214 *
 215 *   int - TRUE if intersections were found.
 216 *
 217 * AUTHOR
 218 *
 219 *   Alexander Enzmann, Dieter Bayer
 220 *
 221 * DESCRIPTION
 222 *
 223 *   Determine ray/superellipsoid intersection.
 224 *
 225 * CHANGES
 226 *
 227 *   Oct 1994 : Creation.
 228 *
 229 ******************************************************************************/
 230
 231 static int intersect_superellipsoid(RAY *Ray, SUPERELLIPSOID *Superellipsoid, ISTACK *Depth_Stack)
 232 {
 233   int i, cnt, Found = FALSE;
 234   DBL dists[PLANECOUNT+2];
 235   DBL t, t1, t2, v0, v1, len;
 236   VECTOR P, D, P0, P1, P2, P3;
 237
 238   /* Transform the ray into the superellipsoid space. */
 239
 240   MInvTransPoint(P, Ray->Initial, Superellipsoid->Trans);
 241
 242   MInvTransDirection(D, Ray->Direction, Superellipsoid->Trans);
 243
 244   VLength(len, D);
 245
 246   VInverseScaleEq(D, len);
 247
 248   /* Intersect superellipsoid's bounding box. */
 249
 250   if (!intersect_box(P, D, &t1, &t2))
 251   {
 252     return(FALSE);
 253   }
 254
 255   /* Test if superellipsoid lies 'behind' the ray origin. */
 256
 257   if (t2 < DEPTH_TOLERANCE)
 258   {
 259     return(FALSE);
 260   }
 261
 262   cnt = 0;
 263
 264   if (t1 < DEPTH_TOLERANCE)
 265   {
 266     t1 = DEPTH_TOLERANCE;
 267   }
 268
 269   dists[cnt++] = t1;
 270   dists[cnt++] = t2;
 271
 272   /* Intersect ray with planes cutting superellipsoids in pieces. */
 273
 274   cnt = find_ray_plane_points(P, D, cnt, dists, t1, t2);
 275
 276   if (cnt <= 1)
 277   {
 278     return(FALSE);
 279   }
 280
 281   VEvaluateRay(P0, P, dists[0], D)
 282
 283   v0 = evaluate_superellipsoid(P0, Superellipsoid);
 284
 285   if (fabs(v0) < ZERO_TOLERANCE)
 286   {
 287     if (insert_hit((OBJECT *)Superellipsoid, Ray, dists[0] / len, Depth_Stack))
 288     {
 289       if (Superellipsoid->Type & IS_CHILD_OBJECT)
 290       {
 291         Found = TRUE;
 292       }
 293       else
 294       {
 295         return(TRUE);
 296       }
 297     }
 298   }
 299
 300   for (i = 1; i < cnt; i++)
 301   {
 302     VEvaluateRay(P1, P, dists[i], D)
 303
 304     v1 = evaluate_superellipsoid(P1, Superellipsoid);
 305
 306     if (fabs(v1) < ZERO_TOLERANCE)
 307     {
 308       if (insert_hit((OBJECT *)Superellipsoid, Ray, dists[i] / len, Depth_Stack))
 309       {
 310         if (Superellipsoid->Type & IS_CHILD_OBJECT)
 311         {
 312           Found = TRUE;
 313         }
 314         else
 315         {
 316           return(TRUE);
 317         }
 318       }
 319     }
 320     else
 321     {
 322       if (v0 * v1 < 0.0)
 323       {
 324         /* Opposite signs, there must be a root between */
 325
 326         solve_hit1(Superellipsoid, v0, P0, v1, P1, P2);
 327
 328         VSub(P3, P2, P);
 329
 330         VLength(t, P3);
 331
 332         if (insert_hit((OBJECT *)Superellipsoid, Ray, t / len, Depth_Stack))
 333         {
 334           if (Superellipsoid->Type & IS_CHILD_OBJECT)
 335           {
 336             Found = TRUE;
 337           }
 338           else
 339           {
 340             return(TRUE);
 341           }
 342         }
 343       }
 344       else
 345       {
 346         /*
 347          * Although there was no sign change, we may actually be approaching
 348          * the surface. In this case, we are being fooled by the shape of the
 349          * surface into thinking there isn't a root between sample points.
 350          */
 351
 352         if (check_hit2(Superellipsoid, P, D, dists[i-1], P0, v0, dists[i], &t, P2))
 353         {
 354           if (insert_hit((OBJECT *)Superellipsoid, Ray, t / len, Depth_Stack))
 355           {
 356             if (Superellipsoid->Type & IS_CHILD_OBJECT)
 357             {
 358               Found = TRUE;
 359             }
 360             else
 361             {
 362               return(TRUE);
 363             }
 364           }
 365           else
 366           {
 367             break;
 368           }
 369         }
 370       }
 371     }
 372
 373     v0 = v1;
 374
 375     Assign_Vector(P0, P1);
 376   }
 377
 378   return(Found);
 379 }
 380
 381
 382
 383 /*****************************************************************************
 384 *
 385 * FUNCTION
 386 *
 387 *   Inside_Superellipsoid
 388 *
 389 * INPUT
 390 *
 391 *   IPoint - Intersection point
 392 *   Object - Object
 393 *
 394 * OUTPUT
 395 *
 396 * RETURNS
 397 *
 398 *   int - TRUE if inside
 399 *
 400 * AUTHOR
 401 *
 402 *   Dieter Bayer
 403 *
 404 * DESCRIPTION
 405 *
 406 *   Test if a point lies inside the superellipsoid.
 407 *
 408 * CHANGES
 409 *
 410 *   Oct 1994 : Creation.
 411 *
 412 ******************************************************************************/
 413
 414 static int Inside_Superellipsoid(VECTOR IPoint, OBJECT *Object)
 415 {
 416   DBL val;
 417   VECTOR P;
 418   SUPERELLIPSOID *Superellipsoid = (SUPERELLIPSOID *)Object;
 419
 420   /* Transform the point into the superellipsoid space. */
 421
 422   MInvTransPoint(P, IPoint, Superellipsoid->Trans);
 423
 424   val = evaluate_superellipsoid(P, Superellipsoid);
 425
 426   if (val < EPSILON)
 427   {
 428     return(!Test_Flag(Superellipsoid, INVERTED_FLAG));
 429   }
 430   else
 431   {
 432     return(Test_Flag(Superellipsoid, INVERTED_FLAG));
 433   }
 434 }
 435
 436
 437
 438 /*****************************************************************************
 439 *
 440 * FUNCTION
 441 *
 442 *   Superellipsoid_Normal
 443 *
 444 * INPUT
 445 *
 446 *   Result - Normal vector
 447 *   Object - Object
 448 *   Inter  - Intersection found
 449 *
 450 * OUTPUT
 451 *
 452 *   Result
 453 *
 454 * RETURNS
 455 *
 456 * AUTHOR
 457 *
 458 *   Dieter Bayer
 459 *
 460 * DESCRIPTION
 461 *
 462 *   Calculate the normal of the superellipsoid in a given point.
 463 *
 464 * CHANGES
 465 *
 466 *   Oct 1994 : Creation.
 467 *
 468 ******************************************************************************/
 469
 470 static void Superellipsoid_Normal(VECTOR Result, OBJECT *Object, INTERSECTION *Inter)
 471 {
 472   DBL k, x, y, z;
 473   VECTOR P, N, E;
 474   SUPERELLIPSOID *Superellipsoid = (SUPERELLIPSOID *)Object;
 475
 476   /* Transform the point into the superellipsoid space. */
 477
 478   MInvTransPoint(P, Inter->IPoint, Superellipsoid->Trans);
 479
 480   Assign_Vector(E, Superellipsoid->Power);
 481
 482   x = fabs(P[X]);
 483   y = fabs(P[Y]);
 484   z = fabs(P[Z]);
 485
 486   k = power(power(x, E[X]) + power(y, E[X]), E[Y] - 1.0);
 487
 488   N[X] = k * SGNX(P[X]) * power(x, E[X] - 1.0);
 489   N[Y] = k * SGNX(P[Y]) * power(y, E[X] - 1.0);
 490   N[Z] =     SGNX(P[Z]) * power(z, E[Z] - 1.0);
 491
 492   /* Transform the normalt out of the superellipsoid space. */
 493
 494   MTransNormal(Result, N, Superellipsoid->Trans);
 495
 496   VNormalize(Result, Result);
 497 }
 498
 499
 500
 501 /*****************************************************************************
 502 *
 503 * FUNCTION
 504 *
 505 *   Translate_Superellipsoid
 506 *
 507 * INPUT
 508 *
 509 *   Object - Object
 510 *   Vector - Translation vector
 511 *
 512 * OUTPUT
 513 *
 514 *   Object
 515 *
 516 * RETURNS
 517 *
 518 * AUTHOR
 519 *
 520 *   Dieter Bayer
 521 *
 522 * DESCRIPTION
 523 *
 524 *   Translate a superellipsoid.
 525 *
 526 * CHANGES
 527 *
 528 *   Oct 1994 : Creation.
 529 *
 530 ******************************************************************************/
 531
 532 static void Translate_Superellipsoid(OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
 533 {
 534   Transform_Superellipsoid(Object, Trans);
 535 }
 536
 537
 538
 539 /*****************************************************************************
 540 *
 541 * FUNCTION
 542 *
 543 *   Rotate_Superellipsoid
 544 *
 545 * INPUT
 546 *
 547 *   Object - Object
 548 *   Vector - Rotation vector
 549 *
 550 * OUTPUT
 551 *
 552 *   Object
 553 *
 554 * RETURNS
 555 *
 556 * AUTHOR
 557 *
 558 *   Dieter Bayer
 559 *
 560 * DESCRIPTION
 561 *
 562 *   Rotate a superellipsoid.
 563 *
 564 * CHANGES
 565 *
 566 *   Oct 1994 : Creation.
 567 *
 568 ******************************************************************************/
 569
 570 static void Rotate_Superellipsoid(OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
 571 {
 572   Transform_Superellipsoid(Object, Trans);
 573 }
 574
 575
 576
 577 /*****************************************************************************
 578 *
 579 * FUNCTION
 580 *
 581 *   Scale_Superellipsoid
 582 *
 583 * INPUT
 584 *
 585 *   Object - Object
 586 *   Vector - Scaling vector
 587 *
 588 * OUTPUT
 589 *
 590 *   Object
 591 *
 592 * RETURNS
 593 *
 594 * AUTHOR
 595 *
 596 *   Dieter Bayer
 597 *
 598 * DESCRIPTION
 599 *
 600 *   Scale a superellipsoid.
 601 *
 602 * CHANGES
 603 *
 604 *   Oct 1994 : Creation.
 605 *
 606 ******************************************************************************/
 607
 608 static void Scale_Superellipsoid(OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
 609 {
 610   Transform_Superellipsoid(Object, Trans);
 611 }
 612
 613
 614
 615 /*****************************************************************************
 616 *
 617 * FUNCTION
 618 *
 619 *   Transform_Superellipsoid
 620 *
 621 * INPUT
 622 *
 623 *   Object - Object
 624 *   Trans  - Transformation to apply
 625 *
 626 * OUTPUT
 627 *
 628 *   Object
 629 *
 630 * RETURNS
 631 *
 632 * AUTHOR
 633 *
 634 *   Dieter Bayer
 635 *
 636 * DESCRIPTION
 637 *
 638 *   Transform a superellipsoid and recalculate its bounding box.
 639 *
 640 * CHANGES
 641 *
 642 *   Oct 1994 : Creation.
 643 *
 644 ******************************************************************************/
 645
 646 static void Transform_Superellipsoid(OBJECT *Object, TRANSFORM *Trans)
 647 {
 648   Compose_Transforms(((SUPERELLIPSOID *)Object)->Trans, Trans);
 649
 650   Compute_Superellipsoid_BBox((SUPERELLIPSOID *)Object);
 651 }
 652
 653
 654
 655 /*****************************************************************************
 656 *
 657 * FUNCTION
 658 *
 659 *   Invert_Superellipsoid
 660 *
 661 * INPUT
 662 *
 663 *   Object - Object
 664 *
 665 * OUTPUT
 666 *
 667 *   Object
 668 *
 669 * RETURNS
 670 *
 671 * AUTHOR
 672 *
 673 *   Dieter Bayer
 674 *
 675 * DESCRIPTION
 676 *
 677 *   Invert a superellipsoid.
 678 *
 679 * CHANGES
 680 *
 681 *   Oct 1994 : Creation.
 682 *
 683 ******************************************************************************/
 684
 685 static void Invert_Superellipsoid(OBJECT *Object)
 686 {
 687   Invert_Flag(Object, INVERTED_FLAG);
 688 }
 689
 690
 691
 692 /*****************************************************************************
 693 *
 694 * FUNCTION
 695 *
 696 *   Create_Superellipsoid
 697 *
 698 * INPUT
 699 *
 700 * OUTPUT
 701 *
 702 * RETURNS
 703 *
 704 *   SUPERELLIPSOID * - new superellipsoid
 705 *
 706 * AUTHOR
 707 *
 708 *   Dieter Bayer
 709 *
 710 * DESCRIPTION
 711 *
 712 *   Create a new superellipsoid.
 713 *
 714 * CHANGES
 715 *
 716 *   Oct 1994 : Creation.
 717 *
 718 ******************************************************************************/
 719
 720 SUPERELLIPSOID *Create_Superellipsoid()
 721 {
 722   SUPERELLIPSOID *New;
 723
 724   New = (SUPERELLIPSOID *)POV_MALLOC(sizeof(SUPERELLIPSOID), "superellipsoid");
 725
 726   INIT_OBJECT_FIELDS(New,SUPERELLIPSOID_OBJECT,&Superellipsoid_Methods)
 727
 728   New->Trans = Create_Transform();
 729
 730   Make_Vector(New->Power, 2.0, 2.0, 2.0);
 731
 732   return(New);
 733 }
 734
 735
 736
 737 /*****************************************************************************
 738 *
 739 * FUNCTION
 740 *
 741 *   Copy_Superellipsoid
 742 *
 743 * INPUT
 744 *
 745 *   Object - Object
 746 *
 747 * OUTPUT
 748 *
 749 * RETURNS
 750 *
 751 *   void * - New superellipsoid
 752 *
 753 * AUTHOR
 754 *
 755 *   Dieter Bayer
 756 *
 757 * DESCRIPTION
 758 *
 759 *   Copy a superellipsoid structure.
 760 *
 761 * CHANGES
 762 *
 763 *   Oct 1994 : Creation.
 764 *
 765 ******************************************************************************/
 766
 767 static SUPERELLIPSOID *Copy_Superellipsoid(OBJECT *Object)
 768 {
 769   SUPERELLIPSOID *New, *Superellipsoid = (SUPERELLIPSOID *)Object;
 770
 771   New = Create_Superellipsoid();
 772
 773   /* Get rid of the transformation created in Create_Superellipsoid(). */
 774
 775   Destroy_Transform(New->Trans);
 776
 777   /* Copy superellipsoid. */
 778
 779   *New = *Superellipsoid;
 780
 781   New->Trans = Copy_Transform(Superellipsoid->Trans);
 782
 783   return(New);
 784 }
 785
 786
 787
 788 /*****************************************************************************
 789 *
 790 * FUNCTION
 791 *
 792 *   Destroy_Superellipsoid
 793 *
 794 * INPUT
 795 *
 796 *   Object - Object
 797 *
 798 * OUTPUT
 799 *
 800 *   Object
 801 *
 802 * RETURNS
 803 *
 804 * AUTHOR
 805 *
 806 *   Dieter Bayer
 807 *
 808 * DESCRIPTION
 809 *
 810 *   Destroy a superellipsoid.
 811 *
 812 * CHANGES
 813 *
 814 *   Oct 1994 : Creation.
 815 *
 816 ******************************************************************************/
 817
 818 static void Destroy_Superellipsoid(OBJECT *Object)
 819 {
 820   SUPERELLIPSOID *Superellipsoid = (SUPERELLIPSOID *)Object;
 821
 822   Destroy_Transform(Superellipsoid->Trans);
 823
 824   POV_FREE (Object);
 825 }
 826
 827
 828
 829 /*****************************************************************************
 830 *
 831 * FUNCTION
 832 *
 833 *   Compute_Superellipsoid_BBox
 834 *
 835 * INPUT
 836 *
 837 *   Superellipsoid - Superellipsoid
 838 *
 839 * OUTPUT
 840 *
 841 *   Superellipsoid
 842 *
 843 * RETURNS
 844 *
 845 * AUTHOR
 846 *
 847 *   Dieter Bayer
 848 *
 849 * DESCRIPTION
 850 *
 851 *   Calculate the bounding box of a superellipsoid.
 852 *
 853 * CHANGES
 854 *
 855 *   Oct 1994 : Creation.
 856 *
 857 ******************************************************************************/
 858
 859 void Compute_Superellipsoid_BBox(SUPERELLIPSOID *Superellipsoid)
 860 {
 861   Make_BBox(Superellipsoid->BBox, -1.0001, -1.0001, -1.0001, 2.0002, 2.0002, 2.0002);
 862
 863   Recompute_BBox(&Superellipsoid->BBox, Superellipsoid->Trans);
 864 }
 865
 866
 867
 868 /*****************************************************************************
 869 *
 870 * FUNCTION
 871 *
 872 *   intersect_box
 873 *
 874 * INPUT
 875 *
 876 *   P, D       - Ray origin and direction
 877 *   dmin, dmax - Intersection depths
 878 *
 879 * OUTPUT
 880 *
 881 *   dmin, dmax
 882 *
 883 * RETURNS
 884 *
 885 *   int - TRUE, if hit
 886 *
 887 * AUTHOR
 888 *
 889 *   Dieter Bayer
 890 *
 891 * DESCRIPTION
 892 *
 893 *   Intersect a ray with an axis aligned unit box.
 894 *
 895 * CHANGES
 896 *
 897 *   Oct 1994 : Creation.
 898 *
 899 ******************************************************************************/
 900
 901 static int intersect_box(VECTOR P, VECTOR  D, DBL *dmin, DBL  *dmax)
 902 {
 903   DBL tmin = 0.0, tmax = 0.0;
 904
 905   /* Left/right. */
 906
 907   if (fabs(D[X]) > EPSILON)
 908   {
 909     if (D[X] > EPSILON)
 910     {
 911       *dmin = (MIN_VALUE - P[X]) / D[X];
 912
 913       *dmax = (MAX_VALUE - P[X]) / D[X];
 914
 915       if (*dmax < EPSILON) return(FALSE);
 916     }
 917     else
 918     {
 919       *dmax = (MIN_VALUE - P[X]) / D[X];
 920
 921       if (*dmax < EPSILON) return(FALSE);
 922
 923       *dmin = (MAX_VALUE - P[X]) / D[X];
 924     }
 925
 926     if (*dmin > *dmax) return(FALSE);
 927   }
 928   else
 929   {
 930     if ((P[X] < MIN_VALUE) || (P[X] > MAX_VALUE))
 931     {
 932       return(FALSE);
 933     }
 934
 935     *dmin = -BOUND_HUGE;
 936     *dmax =  BOUND_HUGE;
 937   }
 938
 939   /* Top/bottom. */
 940
 941   if (fabs(D[Y]) > EPSILON)
 942   {
 943     if (D[Y] > EPSILON)
 944     {
 945       tmin = (MIN_VALUE - P[Y]) / D[Y];
 946
 947       tmax = (MAX_VALUE - P[Y]) / D[Y];
 948     }
 949     else
 950     {
 951       tmax = (MIN_VALUE - P[Y]) / D[Y];
 952
 953       tmin = (MAX_VALUE - P[Y]) / D[Y];
 954     }
 955
 956     if (tmax < *dmax)
 957     {
 958       if (tmax < EPSILON) return(FALSE);
 959
 960       if (tmin > *dmin)
 961       {
 962         if (tmin > tmax) return(FALSE);
 963
 964         *dmin = tmin;
 965       }
 966       else
 967       {
 968         if (*dmin > tmax) return(FALSE);
 969       }
 970
 971       *dmax = tmax;
 972     }
 973     else
 974     {
 975       if (tmin > *dmin)
 976       {
 977         if (tmin > *dmax) return(FALSE);
 978
 979         *dmin = tmin;
 980       }
 981     }
 982   }
 983   else
 984   {
 985     if ((P[Y] < MIN_VALUE) || (P[Y] > MAX_VALUE))
 986     {
 987       return(FALSE);
 988     }
 989   }
 990
 991   /* Front/back. */
 992
 993   if (fabs(D[Z]) > EPSILON)
 994   {
 995     if (D[Z] > EPSILON)
 996     {
 997       tmin = (MIN_VALUE - P[Z]) / D[Z];
 998
 999       tmax = (MAX_VALUE - P[Z]) / D[Z];
1000     }
1001     else
1002     {
1003       tmax = (MIN_VALUE - P[Z]) / D[Z];
1004
1005       tmin = (MAX_VALUE - P[Z]) / D[Z];
1006     }
1007
1008     if (tmax < *dmax)
1009     {
1010       if (tmax < EPSILON) return(FALSE);
1011
1012       if (tmin > *dmin)
1013       {
1014         if (tmin > tmax) return(FALSE);
1015
1016         *dmin = tmin;
1017       }
1018       else
1019       {
1020         if (*dmin > tmax) return(FALSE);
1021       }
1022
1023       *dmax = tmax;
1024     }
1025     else
1026     {
1027       if (tmin > *dmin)
1028       {
1029         if (tmin > *dmax) return(FALSE);
1030
1031         *dmin = tmin;
1032       }
1033     }
1034   }
1035   else
1036   {
1037     if ((P[Z] < MIN_VALUE) || (P[Z] > MAX_VALUE))
1038     {
1039       return(FALSE);
1040     }
1041   }
1042
1043   *dmin = tmin;
1044   *dmax = tmax;
1045
1046   return(TRUE);
1047 }
1048
1049
1050
1051 /*****************************************************************************
1052 *
1053 * FUNCTION
1054 *
1055 *   evaluate_superellipsoid
1056 *
1057 * INPUT
1058 *
1059 *   P          - Point to evaluate
1060 *   Superellipsoid - Pointer to superellipsoid
1061 *
1062 * OUTPUT
1063 *
1064 * RETURNS
1065 *
1066 *   DBL
1067 *
1068 * AUTHOR
1069 *
1070 *   Dieter Bayer
1071 *
1072 * DESCRIPTION
1073 *
1074 *   Get superellipsoid value in the given point.
1075 *
1076 * CHANGES
1077 *
1078 *   Oct 1994 : Creation.
1079 *
1080 ******************************************************************************/
1081
1082 static DBL evaluate_superellipsoid(VECTOR P, SUPERELLIPSOID *Superellipsoid)
1083 {
1084   VECTOR V1;
1085
1086   V1[X] = power(fabs(P[X]), Superellipsoid->Power[X]);
1087   V1[Y] = power(fabs(P[Y]), Superellipsoid->Power[X]);
1088   V1[Z] = power(fabs(P[Z]), Superellipsoid->Power[Z]);
1089
1090   return(power(V1[X] + V1[Y], Superellipsoid->Power[Y]) + V1[Z] - 1.0);
1091 }
1092
1093
1094
1095 /*****************************************************************************
1096 *
1097 * FUNCTION
1098 *
1099 *   power
1100 *
1101 * INPUT
1102 *
1103 *   x - Argument
1104 *   e - Power
1105 *
1106 * OUTPUT
1107 *
1108 * RETURNS
1109 *
1110 *   DBL
1111 *
1112 * AUTHOR
1113 *
1114 *   Dieter Bayer
1115 *
1116 * DESCRIPTION
1117 *
1118 *   Raise x to the power of e.
1119 *
1120 * CHANGES
1121 *
1122 *   Oct 1994 : Creation.
1123 *
1124 ******************************************************************************/
1125
1126 static DBL power(DBL x, DBL  e)
1127 {
1128   register int i;
1129   register DBL b;
1130
1131   b = x;
1132
1133   i = (int)e;
1134
1135   /* Test if we have an integer power. */
1136
1137   if (e == (DBL)i)
1138   {
1139     switch (i)
1140     {
1141       case 0: return(1.0);
1142
1143       case 1: return(b);
1144
1145       case 2: return(Sqr(b));
1146
1147       case 3: return(Sqr(b) * b);
1148
1149       case 4: b *= b; return(Sqr(b));
1150
1151       case 5: b *= b; return(Sqr(b) * x);
1152
1153       case 6: b *= b; return(Sqr(b) * b);
1154
1155       default: return(pow(x, e));
1156     }
1157   }
1158   else
1159   {
1160     return(pow(x, e));
1161   }
1162 }
1163
1164
1165
1166 /*****************************************************************************
1167 *
1168 * FUNCTION
1169 *
1170 *   insert_hit
1171 *
1172 * INPUT
1173 *
1174 *   Object      - Object
1175 *   Ray         - Ray
1176 *   Depth       - Intersection depth
1177 *   Depth_Stack - Intersection stack
1178 *
1179 * OUTPUT
1180 *
1181 *   Depth_Stack
1182 *
1183 * RETURNS
1184 *
1185 *   int - TRUE, if intersection is valid
1186 *
1187 * AUTHOR
1188 *
1189 *   Dieter Bayer
1190 *
1191 * DESCRIPTION
1192 *
1193 *   Push an intersection onto the depth stack if it is valid.
1194 *
1195 * CHANGES
1196 *
1197 *   Nov 1994 : Creation.
1198 *
1199 ******************************************************************************/
1200
1201 static int insert_hit(OBJECT *Object, RAY *Ray, DBL Depth, ISTACK *Depth_Stack)
1202 {
1203   VECTOR IPoint;
1204
1205   if ((Depth > DEPTH_TOLERANCE) && (Depth < Max_Distance))
1206   {
1207     VEvaluateRay(IPoint, Ray->Initial, Depth, Ray->Direction);
1208
1209     if (Point_In_Clip(IPoint, Object->Clip))
1210     {
1211       push_entry(Depth, IPoint, Object, Depth_Stack);
1212
1213       return(TRUE);
1214     }
1215   }
1216
1217   return(FALSE);
1218 }
1219
1220
1221
1222 /*****************************************************************************
1223 *
1224 * FUNCTION
1225 *
1226 *   compdists
1227 *
1228 * INPUT
1229 *
1230 * OUTPUT
1231 *
1232 * RETURNS
1233 *
1234 * AUTHOR
1235 *
1236 *   Alexander Enzmann
1237 *
1238 * DESCRIPTION
1239 *
1240 *   Compare two slabs.
1241 *
1242 * CHANGES
1243 *
1244 *   Nov 1994 : Creation.
1245 *
1246 ******************************************************************************/
1247
1248 static int compdists(CONST void *in_a, CONST void *in_b)
1249 {
1250   DBL a, b;
1251
1252   a = *((DBL *)in_a);
1253   b = *((DBL *)in_b);
1254
1255   if (a < b)
1256   {
1257     return(-1);
1258   }
1259
1260   if (a == b)
1261   {
1262     return(0);
1263   }
1264   else
1265   {
1266     return(1);
1267   }
1268 }
1269
1270
1271
1272 /*****************************************************************************
1273 *
1274 * FUNCTION
1275 *
1276 *   find_ray_plane_points
1277 *
1278 * INPUT
1279 *
1280 * OUTPUT
1281 *
1282 * RETURNS
1283 *
1284 * AUTHOR
1285 *
1286 *   Alexander Enzmann
1287 *
1288 * DESCRIPTION
1289 *
1290 *   Find all the places where the ray intersects the set of
1291 *   subdividing planes through the superquadric.  Return the
1292 *   number of valid hits (within the bounding box).
1293 *
1294 * CHANGES
1295 *
1296 *   Nov 1994 : Creation.
1297 *
1298 ******************************************************************************/
1299
1300 static int find_ray_plane_points(VECTOR P, VECTOR  D, int cnt, DBL *dists, DBL  mindist, DBL  maxdist)
1301 {
1302   int i;
1303   DBL t, d;
1304
1305   /* Since min and max dist are the distance to two of the bounding planes
1306      we are considering, there is a high probablity of missing them due to
1307      round off error. Therefore we adjust min and max. */
1308
1309   t = EPSILON * (maxdist - mindist);
1310
1311   mindist -= t;
1312   maxdist += t;
1313
1314   /* Check the sets of planes that cut apart the superquadric. */
1315
1316   for (i = 0; i < PLANECOUNT; i++)
1317   {
1318     d = (D[0] * planes[i][0] + D[1] * planes[i][1] + D[2] * planes[i][2]);
1319
1320     if (fabs(d) < EPSILON)
1321     {
1322       /* Can't possibly get a hit for this combination of ray and plane. */
1323
1324       continue;
1325     }
1326
1327     t = (planes[i][3] - (P[0] * planes[i][0] + P[1] * planes[i][1] + P[2] * planes[i][2])) / d;
1328
1329     if ((t >= mindist) && (t <= maxdist))
1330     {
1331       dists[cnt++] = t;
1332     }
1333   }
1334
1335   /* Sort the results for further processing. */
1336
1337   QSORT((void *)(dists), (size_t)cnt, sizeof(DBL), compdists);
1338
1339   return(cnt);
1340 }
1341
1342
1343
1344 /*****************************************************************************
1345 *
1346 * FUNCTION
1347 *
1348 *   solve_hit1
1349 *
1350 * INPUT
1351 *
1352 * OUTPUT
1353 *
1354 * RETURNS
1355 *
1356 * AUTHOR
1357 *
1358 *   Alexander Enzmann
1359 *
1360 * DESCRIPTION
1361 *
1362 *   Home in on the root of a superquadric using a combination of
1363 *   secant and bisection methods.  This routine requires that
1364 *   the sign of the function be different at P0 and P1, it will
1365 *   fail drastically if this isn't the case.
1366 *
1367 * CHANGES
1368 *
1369 *   Nov 1994 : Creation.
1370 *
1371 ******************************************************************************/
1372
1373 static void solve_hit1(SUPERELLIPSOID *Super, DBL v0, VECTOR tP0, DBL  v1, VECTOR  tP1, VECTOR  P)
1374 {
1375   int i;
1376   DBL x, v2, v3;
1377   VECTOR P0, P1, P2, P3;
1378
1379   Assign_Vector(P0, tP0);
1380   Assign_Vector(P1, tP1);
1381
1382   /* The sign of v0 and v1 changes between P0 and P1, this
1383      means there is an intersection point in there somewhere. */
1384
1385   for (i = 0; i < MAX_ITERATIONS; i++)
1386   {
1387     if (fabs(v0) < ZERO_TOLERANCE)
1388     {
1389       /* Near point is close enough to an intersection - just use it. */
1390
1391       Assign_Vector(P, P0);
1392
1393       break;
1394     }
1395
1396     if (fabs(v1) < ZERO_TOLERANCE)
1397     {
1398       /* Far point is close enough to an intersection. */
1399
1400       Assign_Vector(P, P1);
1401
1402       break;
1403     }
1404
1405     /* Look at the chord connecting P0 and P1. */
1406
1407     /* Assume a line between the points. */
1408
1409     x = fabs(v0) / fabs(v1 - v0);
1410
1411     VSub(P2, P1, P0);
1412
1413     VAddScaled(P2, P0, x, P2);
1414
1415     v2 = evaluate_superellipsoid(P2, Super);
1416
1417     /* Look at the midpoint between P0 and P1. */
1418
1419     VSub(P3, P1, P0);
1420
1421     VAddScaled(P3, P0, 0.5, P3);
1422
1423     v3 = evaluate_superellipsoid(P3, Super);
1424
1425     if (v0 * v2 > 0.0)
1426     {
1427       if (v1 * v2 > 0.0)
1428       {
1429         /* This should be impossible, since v0 and v1 were opposite signs,
1430            v2 must be either 0 or opposite in sign to either v0 or v1. */
1431
1432         Error("internal failure in function solve_sq_hit1: %d, %g, %g, %g", i, v0, v1, v2);
1433       }
1434       else
1435       {
1436         if (v0 * v3 > 0.0)
1437         {
1438           if (x < 0.5)
1439           {
1440             Assign_Vector(P0, P3);
1441           }
1442           else
1443           {
1444             Assign_Vector(P0, P2);
1445           }
1446         }
1447         else
1448         {
1449           /* We can move both ends. */
1450
1451           Assign_Vector(P0, P2);
1452
1453           Assign_Vector(P1, P3);
1454         }
1455       }
1456     }
1457     else
1458     {
1459       if (v0 * v3 > 0.0)
1460       {
1461         /* We can move both ends. */
1462
1463         Assign_Vector(P0, P3);
1464
1465         Assign_Vector(P1, P2);
1466       }
1467       else
1468       {
1469         if (x < 0.5)
1470         {
1471           Assign_Vector(P1, P2);
1472         }
1473         else
1474         {
1475           Assign_Vector(P1, P3);
1476         }
1477       }
1478     }
1479   }
1480
1481   if (i == MAX_ITERATIONS)
1482   {
1483     /* The loop never quite closed in on the result - just use the point
1484        closest to zero.  This really shouldn't happen since the max number
1485        of iterations is enough to converge with straight bisection.  */
1486
1487     if (fabs(v0) < fabs(v1))
1488     {
1489       Assign_Vector(P, P0);
1490     }
1491     else
1492     {
1493       Assign_Vector(P, P1);
1494     }
1495   }
1496 }
1497
1498
1499
1500
1501 /*****************************************************************************
1502 *
1503 * FUNCTION
1504 *
1505 *   check_hit2
1506 *
1507 * INPUT
1508 *
1509 * OUTPUT
1510 *
1511 * RETURNS
1512 *
1513 * AUTHOR
1514 *
1515 *   Alexander Enzmann
1516 *
1517 * DESCRIPTION
1518 *
1519 *   Try to find the root of a superquadric using Newtons method.
1520 *
1521 * CHANGES
1522 *
1523 *   Nov 1994 : Creation.
1524 *
1525 ******************************************************************************/
1526
1527 static int check_hit2(SUPERELLIPSOID *Super, VECTOR P, VECTOR  D, DBL t0, VECTOR  P0, DBL  v0, DBL  t1, DBL  *t, VECTOR  Q)
1528 {
1529   int i;
1530   DBL dt0, dt1, v1, deltat, maxdelta;
1531   VECTOR P1;
1532
1533   dt0 = t0;
1534   dt1 = t0 + 0.0001 * (t1 - t0);
1535
1536   maxdelta = t1 - t0;
1537
1538   for (i = 0; (dt0 < t1) && (i < MAX_ITERATIONS); i++)
1539   {
1540     VEvaluateRay(P1, P, dt1, D)
1541
1542     v1 = evaluate_superellipsoid(P1, Super);
1543
1544     if (v0 * v1 < 0)
1545     {
1546       /* Found a crossing point, go back and use normal root solving. */
1547
1548       solve_hit1(Super, v0, P0, v1, P1, Q);
1549
1550       VSub(P0, Q, P);
1551
1552       VLength(*t, P0);
1553
1554       return(TRUE);
1555     }
1556     else
1557     {
1558       if (fabs(v1) < ZERO_TOLERANCE)
1559       {
1560          VEvaluateRay(Q, P, dt1, D)
1561
1562          *t = dt1;
1563
1564          return(TRUE);
1565       }
1566       else
1567       {
1568         if (((v0 > 0.0) && (v1 > v0)) || ((v0 < 0.0) && (v1 < v0)))
1569         {
1570           /* We definitely failed */
1571
1572           break;
1573         }
1574         else
1575         {
1576           if (v1 == v0)
1577           {
1578             break;
1579           }
1580           else
1581           {
1582             deltat = v1 * (dt1 - dt0) / (v1 - v0);
1583           }
1584         }
1585       }
1586     }
1587
1588     if (fabs(deltat) > maxdelta)
1589     {
1590        break;
1591     }
1592
1593     v0 = v1;
1594
1595     dt0 = dt1;
1596
1597     dt1 -= deltat;
1598   }
1599
1600   return(FALSE);
1601 }